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A TWO-DIMENSIONAL DIGITAL COMPUTER PROGRAM FOR CALCULATION 


OP OPTIMUM TRAJECTORIES PROM LAUNCH TO INJECTION 

SUMMARY 


A two-dimensional digital computer program has been developed for 
computing optimum trajectories from launch to injection An optimum 
trajectory is the trajectory that allows the maximum payload to be 
placed into the desired orbit Assuming a fixed thrust, the optimum 
trajectory is a minimum flight time trajectory The program equations 
of motion describe the flight of a point-mass multistage vehicle over 
a spherical planet with the effect of planet rotation taken into account 

This program was written to replace two existing trajectory pro- 
grams an atmospheric preset angle -of -attack trajectory program and a 
vacuum calculus of variations trajectory program Preliminary results 
indicate that use of the new program will effect at least a 50 percent 
reduction m the total time required to complete a launch vehicle per- 
formance study using the two existing programs 


INTRODUCTION 


Most preliminary performance trajectory studies necessitate the 
use of large storage, high speed, digital computers Due to the high 
rental cost of these computers, it is highly beneficial to use the most 
efficient programs available Consequently, the Performance Analysis 
Section of the Advanced Spacecraft Technology Division is conducting a 
continuing search for new and improved methods of calculation 

One area in which this continuing study is yielding significant 
progress is that concerned with the calculation of optimum trajectories 
In the past, optimum trajectory calculations have required the use of 
two computer programs* an atmospheric preset angle-of -attack program 
and a vacuum calculus of variations program The former program was 
used to calculate the atmospheric segment of the trajectory as well as 
the initial conditions for the calculus of variations program Using 
these initial conditions, the latter program was subsequently used to 
calculate an optimum trajectory to a desired flight path angle and al- 
titude As used here, an optimum trajectory is that trajectory which 
yields the maximum payload at injection 

The above process for calculating optimum trajectories was time 
consuming in terms of both man hours and machine time Using this pro- 
cess, a minimum of five preset angle-of -attack trajectories were required 
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necessitating 10 separate machine runs with considerable intermediate 
manual data preparation At least 8l calculus of "variations trajec- 
tories were required, or, under less favorable conditions, an average 
of 137 computed trajectories per study 

The program described in this paper was written to replace the 
two programs described above and represents a significant improvement 
over the above technique Using the present program, the average num- 
ber of preset angle-of -attack trajectories computed is seven and the 
average number of calculus of variations trajectories computed is 80, 
all on one computer run This program permits a saving m IBM 709^ 
machine tame of almost 50 percent and a saving in man hours of approxi- 
mately 60 percent 

The present paper presents a detailed description of the new pro- 
gram The appendixes include the equations of motion and auxiliary 
equations, input locations, print schedule, and Fortran listing of the 
program 


SYMBOIS 


Symbol 

Fortran 

Print 

AAA 

AAA 

AM 

A 

AREA 


Ae 

AE 

AEX 

h 

ALT 

ALT 

a 

°D 

CD 

CD 

c* 

CL 

CL 

a 

D 

DRAG 

DRAG 

F 

F 

THRUST 


Definition 
Integration constant 
Cross sectional area of the first 

p 

stage (fu ) 

Exhaust area of the first 
stage engines (ft ) 

Altitude (ft). 

Semi -major axis (ft) 

Axial drag coefficient 

Normal lift coefficient /radian 

Drag (lbs). 

Vehicle thrust (lbs). 
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Symbol 

Fortran 

Print 

Definition 

S 

GRAV 

GRAV 

Vehicle acceleration due to the 

2 

gravity of the planet (ft /sec' ) 

s oe 

GRAVO 


Factor for converting weight to 
mass (ft /sec ) 

I 

sp 

SXP 

ISP 

Specific impulse (see) 

L 

ALIFT 

LIFT 

Lift (lbs). 

m 

SAM 

MASS 

Mass (slugs). 

maeh 

MACH 

MACH 

Mach number 

P 

PRES 

PRES 

Atmospheric pressure (lbs/ft ) 

<2. 

QD 

Q 

Dynamic pressure (lbs/ft 2 ) 

R 

R 

RER 

Distance from center of planet to 
center of vehicle (ft) 

t 

T 

TIME 

Time (sec) 

V 

V 

VEL 

Vehicle velocity (ft/sec) 

w 

WGT 

WLBS 

Weight (lbs) 

X 

X 

XXX 

Ground range 

ct 

ALPHA 

ALPHA 

Angle-of -attack (deg). 

0 

TATER 

THETA 

Angle between local vertical and 
velocity vector (deg) 

X 

ALAM 

LAMB 

Lagrange Multiplier 

P 

RHO 

RHO 

Atmospheric density (lbs/ft^) 

© 

PHI 

PHI 

Range angle (deg) 

i* 

PHIPR 


Latitude (deg ) 

to 

OMEGA 


Rotation velocity of planet 
(rad/sec ) . 
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Subscripts 

DES 

M 

o 

p 

S 

SF 

SL 


Desired 

Maximum 

Initial or planet surface conditions 

Predicted 

Speed of sound 

Space fixed 

Sea level 

PROGRAM DESCRIPTION 


Atmospheric Preset Angle-of -Attack Computations 

The preset angle-of-attack section of the program is designed to 
predict the point mass trajectory through an atmosphere over a non- 
rotating celestial sphere The equations of motion are written assum- 
ing two degrees of freedom with motion m the pitch plane only. These 
equations together with the necessary auxiliary equations are presented 
in appendix A 

The preset angle-of-attack section simulates the vertical flight 
of the vehicle until some preselected time when the vehicle will have 
obtained sufficient altitude to clear all ground obstructions At this 
time a small angle-of-attack is ramped m and held at a constant value 
until it is ramped out at some preselected time before the high dynamic 
pressure region A no-lift (zero angle-of -attack, gravity tilt) trajec- 
tory is then flown through the high dynamic pressure region until cut- 
off time of the first stage 

Capability of instantaneous changes in flow rate, thrust, and ex- 
haust area are incorporated into this section of the program to simulate 
engine cut-off, coast period, engine failure, and other sudden shifts 
that can occur during flight 

The effect of planet rotation on the velocity of the ascending ve- 
hicle is accounted for by a space-fixing equation introduced after cut- 
off of the first stage This rotational effect is estimated by vecto- 
nally adding the rotational velocity vector of the planet to the vehicle 
velocity vector This calculation produces a change m both the magni- 
tude and direction of the vehicle velocity vector 
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The magnitude of the angle-of -attack introduced during the first 
stage flight is the only free variable of the first stage and is con- 
trolled by a maximization routine 

Vacuum Calculus of Variations Computations 

The vacuum calculus of variations section of the program is a two 
degree-of -freedom trajectory simulation which assumes flight in vacuo 
over a non-rotating celestial sphere The equations of motion for this 
section are identical to those for the preset angle-of -attack section 
with the lift and drag terms omitted The Euler -Lagrange equations re- 
sulting from the variational calculus theory, are written m such a man- 
ner as to produce the optimum instantaneous steering angle throughout 
the trajectory This section of the program is employed at the beginning 
of second stage burning at which tame the vehicle has acquired sufficient 
altitude that the effects of atmospheric drag and lift can be considered 
negligible 

After ignition of the second stage engines a preselected amount of 
burn time is allotted for the required stabilization of the vehicle At 
the end of this time, the launch escape system is jettisoned. The re- 
sulting weight loss causes a discontinuity in the state variable The 
Lagrange multipliers are held constant across this discontinuity to in- 
sure an optimum trajectory This can be done since the problem is one 
of minimum time There are no other discontinuities m this section of 
the program 

The Euler -Lagrange equations are concerned with three variables 
The first of these equations results in the free variable A^ which 

can be used to control the range When A^ is set equal to zero, this 

section of the program will compute the range that will allow the maxi- 
mum payload to be placed in the desired orbit. 

The remaining two variables are the angle-of -attack, <x Q , and the 
time rate of change of angle-of -attack, d^, at second stage ignition 

The Euler -Lagrange equations are manipulated so these two variables may 
be satisfied by either a maximization or an isolation scheme. In this 
program, <Xq and o,q are controlled by an isolation routine and no con- 
trol is exercised over the variable A^ 

The calculus of variations trajectory computations are terminated 
when the cut-off equation is satisfied 

The integration throughout the program is performed by the SHARE 
subroutine RW-IKT 



Isolation Bcmtme 


The isolation routine uses the variables o,q, cx,q, and the altitude 

and flight path angle at the tune of satisfaction of the cut-off equation 
to insure injection into the desired orbit For example, if the desired 
orbit is a 100-nautical mile circular orbit, the cut-off equation will be 
satisfied at local circular velocity The initial tune rate of change 
of angle-of -attack will be varied until the desired flight path angle 
(90 deg) is reached, and the initial angle-of -attack will be varied un- 
til the desired altitude (100 nautical miles) is obtained 

This routine is designed to isolate on the desired end conditions 
m a minimum number of trajectories However, m cases where the ini- 
tial conditions are unfavorable, the routine will use as many as five 
points m order to insure an isolation In all cases, the equations of 
isolation are the same with the following exceptions for the isolation 
of flight path angle, the equations of isolation over a Q use all tra- 
jectories, for the isolation of altitude, the equations of isolation over 
o,q use only those trajectories m which the desired path angle is al- 
ready obtained 

An initial guess for and cc-q and the magnitude of their step 

sizes must be fed into the computer as input These guesses must not 
be equal to zero 

The first trajectory is computed using the input values for the 
initial <Xq and A second trajectory is then computed, using for 

its initial conditions the input a Q and the input a Q summed to the 
step to be taken on <x Q That is, a Q is the same for "both trajecto- 
ries, while o-q is changed by a small amount for the second trajectory 
The direction of the step taken on <%q is fixed by the program so that 

the resulting path angle changes in the direction of the desired path 
angle 

After the first two trajectories are run, the two-point mter- 
polation/extrapolation equation 


a 3 = a 2 * ( 0 DES “ 0 2) (e 2 - 0 1 ) ^ 

is used to arrive at the desired a.^ used to compute the third trajec- 
tory The <Xq' s from these three trajectories are then used to construct 
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a conic type curve which closely duplicates the actual plot of <Xq ver- 
sus the cut-off path angle An accurate three-point fit is obtained 
using the conic equation 


-B ± 
a = 


where 


A 


1 H 
CD 

tf 1 

OJ 

d H 

1 


' 1~ 



p 



B 


OJ 

CD 

1 


1 

C 


K"v 1 
CD 

OJ 

J 


1 


V® 2 - 4A ( Ce DES - 0 
2A 


( 2 ) 


(3) 


Equations (2) and (3) are used to determine the cc^ for the next 

trajectory This type of conic fit works equally as well with four or 
five points For example, for four and five points, equation ( 2 ) becomes 


a = 



(M 


where 


A = A (5) 

BB = BB BES + D (6) 

CC = ^DES + E0 DB5 ' 


1 


(7) 
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For four points, equation (3) becomes 



e 5 2 



(8) 


For five points this equation becomes 



The conic curve fit can use no more than five points; therefore, if 
more than five trajectories are required to obtain the desired value of 
a, the last five computed points are retained 

When the path angle obtained is within tolerance of the desired path 
angle, the trajectory is said to be converged on path angle 

After the routine is converged on path angle, the altitude obtained 
from the converged trajectory is checked against the desired altitude 
A trajectory is then calculated using the o,q from the converged path 

angle m conjunction with the input summed to the step to be taken 

on The direction of the step taken on a Q is fixed by the program 
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so that the resulting altitude of the second converged-on path angle 
trajectory changes m the direction of the desired altitude 

The method of isolating on all subsequent converged-on path angle 
trajectories and converged-on altitude trajectories is the same with the 
exception that the first cx^ and a Q are computed using the knowledge 

gained in converging the previous trajectories Also, the step sizes 
and the direction of the steps on these variables are computed as a func- 
tion of the variations m earlier trajectories 

The information computed from the converged-on altitude trajectory 
is fed into the maximization routine and the isolation routine is re- 
initialized 

Maximization. Routine 

The maximization routine optimizes the complete trajectory using the 
preset angle-of -attack of the first stage trajectory and the cut-off weight 
of the isolated second stage trajectory As used here, an optimum tra- 
jectory is that trajectory which will allow the maximum payload to be 
placed into the desired orbit 

The first isolated second stage trajectory is computed from the end 
point of the first stage trajectory This first stage trajectory is com- 
puted using the input preset angle -of -attack estimate for the first stage 
flight A second first stage trajectory is then computed using the input 
angle -of -attack step added to the preset angle -of -attack estimate The 
preset angle -of -attack step for the third isolated trajectory is taken m 
the direction that will result m an increase in cut-off weight m orbit 
Due to the nature of a plot of preset angle -of -attack as a function of 
cutoff weight m orbit, an extrapolation to the maximum cut-off weight is 
extremely dangerous in that it can result in trajectories that cannot be 
advantageously used by the maximization routine Therefore, additional 
trajectories are computed using a fixed step in the preset angle-of -attack 
until the maximum cut-off weight m orbit has been enclosed 

A knowledge of the range of preset angles-of -attack that can be used 
by the vehicle will result in a shorter computer running time However, 
knowledge of this range is not absolutely necessary if small values are 
input for the preset angle-of -attack and the preset angle -of -attack step 

After the maximum has been enclosed, it is seen that a conic equa- 
tion closely simulates the plot of preset angle -of -attack as a function 
of cut-off weight m orbit The conic equation used m this routine is 



-BW + D 
T> 

2A 


( 10 ) 
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where 


- b ±^\/b 2 - 4ae 


ii 

(11) 

P 

)1 

0 

1 

w ro 

£ 

(12) 

b = E - BD/2A 

(15) 

c — — D 2 /4a - 1 

(14) 


and for a three-point fit. 



For a four-point fit. 




11 


For a live-point fit, 



Equation (10) is used to compute the preset angle -of -attack that will 
result in. the maximum cut-off weight in. orbit. If more than five preset 
angle-of -attack trajectories must be computed in order to insure an opti- 
mum trajectory, the maximization routine will retain and use only the five 
points most likely to produce the maximum 

The cut-off weight m orbit is checked against the maximum cut-off 
weight computed by the maximization routine When these weights are with- 
in tolerance, the program is said bo be m a converged state The pro- 
gram will then print a history of the optimum trajectory variables and 
check to see if additional trajectories are to be computed 

Cut-off Equation 

The conditions which terminate a trajectory calculation program are 
among the limiting factors of the program Hence, the cut-off condition 
should be as versatile as possible Since the computed velocity for a 
conic section is the most adaptable to the computation of a local (as a 
function of altitude) velocity requirement, the conic velocity will be 
employed as a cut-off condition 

The equation for the velocity requirement of all conic sections ex- 
cept the hyperbola can be expressed as 



(18) 
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The velocity requirement for a hyperbola is 



( 19 ) 


One of the most common requirements is that the vehicle pass through 
some point in space m a given direction If the point in space is point 
2(R^, 0^) and the instantaneous position of the vehicle is point 1 



fs R 2 \ 

00 \ 

2(R 2 - R 1 ) R 2 sm 2 e 2 

1 R-, I 

OJ 

OJ 

1 

a 

a 

1 

\ 1 / 

S 2 sin 0 2 - it x sm 0.J 


( 20 ) 


When the required point m space is the apogee of an ellipse, then 
R 2 becomes the radius of the apogee of the ellipse, the flight path angle 

becomes 90°, and the velocity requirement becomes 



/gR 2 \ 

0 O 

2(R 2 - El ) R 2 " 

\ E 1 / 

2 2 2 
R 2 - R 1 sm e 1 

V _L / 


( 21 ) 


The most widely used velocity cut-off requirements are local circular 
and local parabolic escape velocity To compute the local circular ve- 
locity requirement, the quantity m the brackets in equation (21 ) must 
equal one, and to compute the local parabolic escape velocity requirement 
the quantity in the orackets must equal two 


Equation (21) can be written m the following form such that the cal- 
culus of variations section of the program can be terminated on time or 
velocity 



where Z^, Z^, 


Zg, and are input constants 
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In order to terminate the trajectory on some preselected time, the 
following -values are input, and the trajectory time is compared to 
CUTOSF 


Z l = 0 


v° 


Z = 0 
3 


Z 2 = 0 


\ = o 


Zg = desired cutoff time 


e 2 = 1 


For velocity termination of the trajectory, the following values 
are input, and the vehicle velocity is compared to CUTOFF 

For local circular 


"i - 1 


Z, = g R 

3 o o 


V° 


R 2 = 1 


z 2 = ° 


z 6 = 0 


For parabolic escape 


g R 

o o 


z 5 = o r 2 = 1 


z 6 - 0 


g R 

o o 


H ' ^He R 2 ' 1 


For elliptic 


z 6 = ° 


g R 
o o 


~ o Rg - Rg 


o 


z 6 = o 


(27) 
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Areas for Future Revision 

The program in its present form is quite limited and a number of 
improvements and expansions are possible These include 

1 An increase in the number of stages that can be calculated using 
the calculus of variations theory 

2 A provision for weight loss table and sea level thrust table to 
be input as a function of time 

3 Addition of the ability to simulate thrust decay 

li- Modification of the isolation and maximization routine such that 
for escape trajectories one variable will be used for isolation and two 
vanabl es will be us^d for maximization 

5 Addition of the ability to optimize the propellant loading of 
each stage to determine the optimum tank size on predesign stages 

6 Revision of the calculus of variations equations so that the 
atmospheric effects may be taken into account 



APPENDIX A 


PRESET ANGLE-OF -ATTACK EQUATIONS 
OF MOTION AND AUXILIARY EQUATIONS 


The atmospheric preset angle-of -attach equations of motion 

R 

X=”V Sine 

R = V cos 9 


, r P D 

V ss — cos a, - — - g cos 0 

m m 

= — sin a + — + sin 0 

iy mv \v R ) 


- F, 


m = 


SL 


I g 
sp oe 


where 


P = P gL + Ae (P o - F) 


( Ro\ 

55 g o \R ) 


D = Aq (C D cos a + a sin a) 
L = Aq (-Cp sm a + a cos a) 


M = 


W_ 


oe 


also 


mach = 


V_ 

V, 


are 

(Al) 

(A2) 

(A3) 

(A4) 

(A5) 

(An) 

(A7) 

(A8) 

(A9) 

(A10) 

(All) 



r 



(A12) 


X = ® + 6 + a (A13) 

Altitude and range angle are given by 

H = R - Ro (AlU) 

and 

* = ^ (A15) 

respectively 

A vector diagram of the summation of the forces and the coordinate 
system is given m figure 1 on page 17 

The preset angle-of -attach program is given by 



t ^ t ^ t n (Al6) 

o 1 

(A17) 



Time 




L8 


Pressure, density and the velocity of sound are table look-up values 
as a function of altitude The drag and lift coefficients are table look- 
up values as a function of Mach number. 

At first stage cut-off, the rotation of the planet is taken into 
account by means of the equations 


+ 2u>R o V cos ® ' 


sm 9 sin A„ + urR ^ cos^ ® ' 
Z o 


(A21) 


and 


e gF = Arc cos 



cos 



(A22) 
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APPENDIX B 

CALCULUS OF VARIATIONS EQUATIONS OF MOTION, 
EULER-LAGRANGE EQUATIONS, AND AUXILIARY EQUATIONS 


The equations of motion for the -vacuum calculus of variations are 
given below The coordinate system and forces are the same as those m 
figure 1 with the exception that in the calculus of variations there are 


no aerodynamic forces acting on the vehicle 



X = V sin 0 

it 

(Bl) 


R = V cos 0 

(B2) 

V s 

F 

= — cos a - g cos 0 
m 

(B3) 

8-*- 

mv 

sm a + (f “ ^) sin 0 

W 


- F 

m = ~ 

sp ®oe 

(B5) 

where 

V = V SF 

( B6) 


e = e SF 

(B7) 

and 

/ R \ 2 

e - s o (r) 

(BB) 


¥ 

m = 

e ° e 

(B9) 
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The Euler-Lagrange equations result m 


H-o 


- zKz 
\ ^ 


R COS 0 + R lv 


X 3 = - \ cos 9 + ~ 


; _ 

5 ' V 3 


V 



(BIO) 

: - i) sxn 6 

(BU) 

; + (v + 1) sin ® 

(B12) 

\ ($ - 1) “• • 

(B13) 

• sm a | 

(B14) 


and 


0 = 1, sm a cos a 

3 v 


a = arc tan 




(B15) 

(B16) 


Differentiating equation (B15) and substituting equation (B-12) and 
(B-13 ) for A^ and leads to 


S = 


A^B cos a 


and 


where 


A^ = A^B sm a. 


(B18) 


(B19) 


B = 


V 2 sin ( a -fr q) 


P v 2 

V a + g sm 0 + ~ sm a - ~ sin a cos (a + e) 


(B20) 



(B2l) 


l 2 = |B| 

It can be shown from equation (£1.4) that 

Xe- ^ 0 when P ^ 0 

Substituting equations (Bl8) and (B19) into equation (Bl4) and 
solving the inequality results m equation (E2l) 

£y solving equations (Bl8) through (B2l) at the time of initiali- 
zation of the calculus of variations the Iagrange multipliers are ex- 
pressed as a function of the state variables and the control variables 
cc Q and a Q 
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APPENDIX C 

INPUT LOCATIONS AND PRINT SCHEDULE 


The subroutine used by this program to input information into the 
machine is RW-FINP The FINP input subroutine uses D, N, J, F, and <t> m 
a unique order and are defined here instead of m the list of symbols 

Ni = D is the introduction of the table being entered into the 
machine where 1 is the number of the table m the calling sequence 
Ji = the introduction of each entry m the table, where 1 is the loca- 
tion number m the table F signals the end of each table and •$> sig- 
nals the end of each calling sequence 

COMMENT CARD 

N1 = D, (Altitude Table) F 
N2 = D, (Density Table) F 
N3 = D, (Pressure Table) F 

= D, (l /Velocity of Sound Table) F^ 

N1 = D, (mach Table) F 
N2 = D, (C D Table) F 


Atmosphere locations, 
r 250 locations per table 

J 


N5 = D, (C L Table) F$ 


Nl = D 

as X q , Initial Range 

J5 = E 5 R + h , Initial Radius to Vehicle 
o o’ 

Jg = V , Initial Velocity 

J„ = 6 F, Initial Flight Path Angle 
0 (from Vertical) 



This block of input may 
be omitted Built m 
values are* 

x = 0 V = 0 

0 o 

R q = 20898906 9=0 

Input values may be used 
to override these and 
other built in values 
(below) 



25 


N2 = D 


J2 



Step size for 1st stage 


This input may he omitted 
Built in. values are* 


J5 = KAJP = IB or 2B 

integration mode for 1st stage 

( 11 ^ 

j4 = At v , initial step size for 
2nd stage 


J5 - KA.IC = 0, IB, or 2B integration 
mode for 2nd stage 




A t « = 
KAIP = 

At™ 
KAJC = 


1 

2B 

* 5 
0 


(integration modes are* 


0 - Variable step Adams -Moulton 


1 - Fixed step Runge-Kutta 

2 - Fixed step Adams -Moulton) 


a 

ii 

Integration Parameters; see writeup on RW-HW 
No input necessary; built m values are 

in = 

A7j 

A2 - ICT 6 , A5 = A4 = A5 

II 

& 

11 

A7 = 0 

J15 = 

g oe 

• ■„ g — conversion factor 
Mass 

-\ 

No input necessary. 

Jl4 = 

s oh 

Surface gravity of body 


built in values are 

J15 a 

R 

0 

Radius of body 


S oe = 32 1849 

J16 = 

»o 

Sea level pressure 

■ 

> s ob = 32 1849 

JIT = 

(0 

Angular velocity of body 


R q = 20898906 

Ji8 = 

$ 1 

latitude 


p = 2124 2l4 
0 

J19 = 

Az 

Azimuth 


-4- 

U) = 7292U5 x 10 


§' = 28 28° 




J21 = Printout Controls; 1 prints every tra- 
jectory according to J22; 0 prints 
only last trajectory 


J22 = nB to print every nth integration 
step when J21 = 1 


STo input necessary, built 
m values are 


J23 ' nB to print every nth integration 
step of final trajectory 




J21 = 0 
J22 = 0 
J23 = 5B 


For both J22 and J23 a 0 causes only 1 st and last points to be printed 


J2b- » ±nB where + or - causes maximization or minimization respectively 
of velocity if n = 1 or weight if n = 2 


J25 = Factor by which J80 is multiplied, 
when a successful step is taken, to 
modify step size 


Built m value 

> 

J25 = 1 


J 28 

J29 

J 30 

J31 

J32 



o 


cut-off time on 1 st stage 
initial weight on 1 st stage 
initial flow rate of 1 st stage 
initial sea level thrust of 1 st stage 
initial exhaust area of 1 st stage 


J33 



cross sectional area for 1 st stage 


J 3 ^ » t , initial time, built m equal to 0 


J35 = 

J38 = t k 


Ramp function times for 
tilt program* 


* Any time not mpux -^ill cause the rest of that particular block to 
be ignored 
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J39 

Jh2 

F^3 

J46 


= t. 


FI 


= t. 


Fk 


F„ 


= F, 


Times during the 1st stage when changes m thrust and 
exhaust area may he specified* 


New values for sea level thrust to he brought in at 
times specified m J39 - J4-2 respectively 


j4T = A a 


New values for exhaust area to he brought in at 
times specified m J39 - j4-2 respectively 


J31 

J5T 


1 

\ Times during the 1st stage when changes in flow 
| rate may he specified * 


J58 

Jo4 

J66 


1 

y New values for flow rate to he brought m at times 
f specified m J51 - J57 respectively 


nB, 2nd stage cut-off control; n = 1 cuts off on velocity, 
n = 2 cuts off on time » J 67 


J 67 



cut-off time of 2nd stage if 


J 66 = 2B 


J68 = TOL, cut-off tolerance 


J69 

J 70 

J 71 

J 72 

J73 

J74 


Z1 

Z2 

Z3 

Z4 

Z5 

z6 


- 0 or 1 
= 0 or 1 

- 8 ob E o 2 
= 0 or 1 

= 0 or 1 

= 0 or 1 




I For velocity cut-off 


J75 - Rg = 1 or apogee radius J 

"Any time not input will cause the rest of that particular block 
to be ignored 



2D 


CUT-OFF EQUATION 


V , = Zl* 
cut 


+ Z5 


(l + 22)* Z3* (l + 2&* R/R 2 )/(l + z4*R 2 /R 2 2 * 3 in 2 e) 


1 

2 


+ z6 


J 76 = a, , initial 1st stage a; 0 if not input 


J77 = a, , maximum a 1st stage tilt program 

lUclX 

J 78 = a Q , initial a of 2nd stage 


J 79 = <x q , initial co of 2 nd stage 


Independent varieties of 
search, initial guesses 
(^ 0 ) must De made 


J 


J80 = to. 


~\ 


max 


j 8 i = Aa S Initial step size or increment for variables to be used 
0 'in search procedure 


J82 = ha. 


O J 


J84 = W initial weight of 2nd stage 

J 85 - IspW, specific impulse of 2nd stage 
J86 = thrust of 2nd stage 

J130 = e DES , desired final flight path angle (from vertical) 
J131 = h^gg, desired final altitude 
J132 = T0L1, tolerance on 0 DEg 
J133 = T0L2, tolerance of h^g 

J13^ = T0L3, tolerance on variable to be maximized 
(velocity or weight) 1$ 
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Print schedule for the preset angle-of -attack computations 


XXX 

XDOT 

TIME 

ALPHA 

ALDOT 

ERR 

RDOT 

ALT 

WLBS 

WDOT 

VEL 

TOOT 

THRUST 

MASS 

MDOT 

THETA 

THDOT 

GRAV 

PHI 

CHI 

RHO 

PRES 

VS 

DRAG 

LIFT 

Q 

AEX 

MACH 

CD 

CL 


Print schedule for the calculus of variations computations 


XXX 

XDOT 

TIME 

ALPHA 

ALDOT 

RRR 

RDOT 

ALT 

WLBS 

WDOT 

VEL 

TOOT 

THRUST 

MASS 

MDOT 

THETA 

THDOT 

GRAV 

PHI 

CHI 

ISP 

LAMES 

LAMB? 

LAMB4 

LAMB5 


DLAM4 


AAA 


DLAM2 


DLAM5 


DLAM5 
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APPENDIX D 


FORTRAN PROGRAM 



29 


SIBFTC MAIN 

CMAIN CONTROL PR06RAM INCLUDING INPUT AND SEARCH PROCEDURE 
DIMENSION TOP ( 12 ) » 0(150)* Tl(150) 

DIMENSION T ( 150 ) »C( 150 ) , ALTA ( 250 > , RHOTA ( 250 ) ,PRETA( 250 ) ,FVSTA( 230 ) 
1»FMTA(30} »CDTA(30) ,CLTA(30) 

DIMENSION C2 ( 150 ) * T2 ( 150 ) 

COMMON T * C 

COMMON ALTA , RHOTA, PRETA ,FVST A , FMTA *CDTA ,CLTA 
COMMON Cl , T1 

COMMON TEST 1 * TEST2 
COMMON C2,T 2 

EQUIVALENCE I TIME* T <2)1 , (STEP,T<3) ),(X,T(4)),(R,T<5)),(V,T(6)), 

1 ( TATER »T 17 ) >,(SAM,T(8) ) » ( ALAM2*T (9 ) ) * ( ALAM3 » T (10) >, (ALAM4,T(11) ) , 

2 < ALAM5 * T ( 12 > ) • ( XDOT * T ( 13 ) ) , ( RDOT , T ( 14 ) ) , (VDOT,T( 15) ),<TATDOT* 

3 TC 16) ) , (SAMDOT *T( 17) ) , (DLAM2,T(18 ) ) , ( DLAM3 , T ( 1 9 ) ) > ( DLAM4 > T ( 20 ) ) , 

4 ( DLAM5 » T 1 2 1 ) ) 

EQUIVALENCE < N ,C ( 1 ) ) , ( STEPP ,C < 2 ) > , C KA1P , C ( 3 ) ) , <STEPC,C<4) ) , (KA1C,C 
1(5) )*(A2*C(6)),(A3,C<7) ) * ( A4 *C ( 8 ) ) , (A5»C(9))»(A6*CllO) ) » ( A? »C ( 11 ) ) 
2 

3» ( GO »C ( 13 > > , ( GRAVO »C( 14 ) ) * ( RO » C ( 15 ) ) * (PRES0.CU6) )» (OMEGA »C< 17) ), < 
4PHIPR,C( 18) ) , (AZ*CI19) ) 

5 

6,(SAVPAP*C(21) ) » ( I PRNT »C < 22 } ) , (NPR,C<23 ) ) »< 10PT ,0 24) ) , (STP,C(25) ) 

7 

8 * ( T I MPC »C ( 28 ) > , <W0*C(29) ) , (WD0.CC30) > . ( FO»C{ 31 ) ) , (AEO»C(32) ) * (AREA 
9 * C ( 33 ) ) , <TIM0*C(34) ) 

EQUIVALENCE < T IM1 , C ( 3 5 ) ) * < TIM2 ,C(36 ) ) , ( T I M3 , C ( 37 ) ) * ( T I M4 ,C < 38 ) I * ( T 
1IMF1 »C(39 ) ) » ( T I MF2 * C ( 40 ) ) , ( T IMF3 . C ( 41 ) ) ♦ ( T1MF4*C<42) ) , < FI ,C (43 > ) , < 
2F2 > C ( 44 ) ) , <F3*C(45> ) » (F4,C(46) ) , CAE1»C(47) ) , (AE2,C148) ) , ( AE3*C(49) 

3) * ( AE4 »C ( 50 ) ) »CT IMW1 » C( 51 ) ) , (T1MW2 ,C( 52 ) 1 » < T I MW3 , C ( 53 ) ),<TIMW4,C(5 
44) ) , ( T IMW5 >C( 55 ) ), (TIMw 6 ,C( 56 > ) , C T IMW7 , C ( 57 ) ) , (UD1»C(58) ) , <WD2,C<5 
59 ) ) » ( WD3 > C ( 60 ) ) » (WD4»C(61) ) ,(WD5.C(62) ) , (WD6*C(63) ) ♦ (#'D7,C(64) ) 

6 

7» ( I CUT, C< 66 ) J , ITIMCC*C<67) ) , (TOL,C<68) ) , ( 21 , C ( 69 ) ) * ( Z2 ,C ( 70 ) ) * ( Z3 , 
8C ( 71 ) ) » ( 24 * C( 72 ) ) * (Z5*C(73) ) > (Z6»C(74) ) > ( R2 > C ( 75 ) ) 

9 

EQUIVALENCE ( ALPHAP,C(76 ) ) , ( ALFMAX , C ( 77 ) ) , ( ALPHAC , C < 78 ) ) , <DALF,C(7 
19) ) , (DLAMX,C(8Q) ) , ( DtLALF , C ( 81 ) ) > ( DELAOT , C < 82 ) ) 

2 

3» < WC»C ( 84 ) ) *(SIP»C(85 ) ) » ( FC *C ( 86 ) ) 

4 

EQUIVALENCE ( ALPHA »C( 90) ) » ( SA»C(91 ) ) , (CA.C (92 ) ) , ( ST»C ( 93 ) ) , ( CT »C ( 9 
14) ) , (SALT»C<95) ) ,( CALT,C<96) ) * <GRAV,C<97> ) » ( ALT >C ( 98 ) ) ,(DRAG»C(99) 
2 ) » ( FVS > C( 100 ) ) »<FMACH,C(101) ) , ( VS, C 1 102 ) ) , ( CD »C { 103 ) ) , (QD ,C ( 104 ) ) , 
3 (RH0»C ( 105 ) ) ><PRES»C( 106 ) ) * C AE , C( 107 ) ) , ( PDIFF , C ( 108 ) ) , { FRAC ,C ( 109 ) 

4) ♦ ( VGT * C ( 110) ) » ( WD0T*C< 111 ) ) » ( F,C( 112 ) ) , ( A, Cl 113) ) * (PHI »C( 114) > , (C 
5HI,C( 115) ) * ( I5AVE1 *Cl 116 ) ) , ( ISAVE2 »C( 117 ) ) » ( JJ , C ( 1 18 ) ) ,IIM0DE,C(11 
69 ) ) > ( U * C( 120 ) ) » (KA1,C( 121) ) » ( CL *C ( 122 ) ) , ( AL I FT , C ( 1 23 ) ) 

7 

8. (TATD,C(130) ) . (ALTD.C! 131 ) ) , ( T0L1 ,C( 132 ) ) » ( T0L2 , C ( 133 ) ) , { T0L3 ,C ( 1 
934) ) 

READ (5*100) TOP 
WRITE (6*100) TOP 
100 FORMAT ( 12A6 ) 

DO 463 1=1,150 
463 C ( I )=0. 

CALL FINP (4, ALTA* RHOTA, PRETA, FVSTA ) 

CALL FINP ( 3*FMTA,CDTA,CLTA ) 

KA1P=2 

KA1C=0 
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STEPP=1. 

STEPC=5. 

A2= 1 . E-6 
T IM0=0» 

ALPHAP=0. 

60=32.1849 

GRAV0=32.1849 

R0=20898906. 

PRES0=2124.214 

0MEGA=.72921l5E-4 

PHIPR=28.28 

A2=90« 

R=20898906. 

X = 0. 

V=0. 

TATER=0. 

SAvPAP=0» 

IPRNT=0 
NPR=5 
STP=1 . 

CAS£2=0 « 

1 CALL FINP (2,T,C) 
N=9 

C0UNT=0. 

Z=57. 295775 
ALFMAX=ALFMAX/Z 
ALPHAC=ALPHAC/Z 
DALF=DALF/Z 
00 2 1=1,150 
C 1 t I )=C( I > 

2 Till >=T< I ) 

T EST 1=0 . 0 
T EST 2=0.0 
T EST 3=0*0 
T EST 4=0 . 

JGOl =1 
JG02=2 
JG03=3 
FAKTOL=0. 

FAKTST=0. 

CASE=0. 

IF (IOPT) 76,76,77 

76 FMNMX=— 1 « 

GO TO 78 

77 FMNMX-1 • 

78 JOPT s 1ABS ( IOPT) 
DADSAV=DELADT/Z 
DALSAV=DELALF/Z 

11 CONTINUE 
AZ=AZ/Z 
PHIPR=PHIPR/Z 
ALPHAP=ALPHAP/Z 
TATER=TATER/Z 
TATD=TATD/Z 
T0L1=T0L1/Z 
D£LALF=DELALF/Z 
DELADT=OELAOT/Z 
DLAMX=DLAMX/Z 

3 CONTINUE 
CALL FCALC 



IF (TEST4}40,16,40 

40 DO 41 1=1*150 
C< I ) = C1 ( I ) 

41 T { I )=T1U J 
WRITE (6*105) 

105 FORMAT <1H1) 

CASE2=1. 

alfmax=alfmax*z 

ALPHAC=ALPHAC*Z 
DALF=DALF*Z 
GO TO 1 

16 CONTINUE 

IF ( TEST1 ) 20*20,21 

20 ALD1 = C2 (79 ) 

T AT1=TATER 

IF ( TEST2 ) 17*17,18 

17 IF ( TEST3 ) 18*19,18 
19 CONTINUE 

IF (CASE2) 80*80*18 

80 CONTINUE 

IF (TATER-TATD) 84*84,72 

84 DELADT=ABS ( DELADT ) 

GO TO 71 

72 DELADT=-ABS< DELADT ) 

GO TO 71 

18 DELADT= ( T ATD-TATER ) *ADI NC/TAT I NC 
IF (DELADT) 71*82*71 

82 DELADT=DADSAV 
GO TO 80 

21 ADINC=C2<79)-ALD1 
T AT I NC=TATER”TAT1 

71 CONTINUE 

COUNT -COUNT+1 « 

IF (COUNT-25*) 300*300,75 
75 WRITE (6*104) 

104 FORMAT (60H025 TRAJECTORIES WITHOUT 1ST ORDER CONVERGENCE • GUESS A 
1GAI N* ) 

GO TO 40 

300 CALL ISO ( TATER *C2 ( 79 ) * TATD ♦ DELADT , TOL1 * TEST 1 * JGOl ) 

IF (TEST1) 3,24,3 
24 IF ( TEST2 ) 26,26,27 

26 ALF1 = C2 (78 ) 

ALT 1=ALT 
DAL 1 = C2 (79 ) 

IF (TEST3) 28,12,28 
12 IF ( CASE2 ) 81,81,28 

81 CONTINUE 

IF (ALT-ALTD) 73,85,85 

85 DELALF=ABS ( DELALF ) 

GO TO 74 

73 DELALF=-ABS ( DELALF ) 

GO TO 74 

28 DELALF=fALTD-ALT>*AL INC/ALT INC 
IF (DELALF) 74,83,74 

83 DELALF=DALSAV 
GO TO 81 

27 ALINC=C2(78)-ALF1 
ALT INC=ALT-ALT1 
DAL INC=C2 ( 79 ) — DALI 
ALFINC=C2(78)— ALF1 
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7 4 CONTINUE 

WRITE (6,102) 

102 FORMAT ( 2 2HOCON VERGED FIRST ORDER) 

COUNT=0. 

ASAV=C2 ( 7 8 ) 

CALL ISO ( ALT »C2 ( 78 ) * ALTD , DELALF »T0L2 »TEST2 , JG02 ) 

IF ( TEST2 ) 14,95,14 
14 CONTINUE 

IF ( TEST3 ) 30,29,30 

29 FAKDEL=0. 

GO TO 4 

30 CONTINUE 

FAKDEL=(C2(78) -ASAV ) #DAL I NC/ALFI NC 

4 CONTINUE 

CALL ISO (ALT, C2(79) , ALTO, FAKDEL ,FAKTOL , FAKTST , JG03 ) 

GO TO 3 
95 CB177=C177 
C177=C1 (77) 

CB278=C278 
C278 = C2 (78 ) 

CB279=C279 

C279=C2C79> 

GO TO (67,68), JOPT 

67 FNC=V 

GO TO 69 

68 FNC=WGT 

69 CONTINUE 

CALL FMXMN < FNC ,C1 < 77 ) , FMNMX*DLAMX ,STP * TOL3 ,TEST3 , JGOl I 
13 WRITE (6,101) 

101 FORMAT (23HOCONVERGED SECOND ORDER) 

FAKTST=0. 

5 DO 6 1=1,150 
CU )=C1CI ) 

6 T( I ) = T 1 ( I ) 

91 IF ( TEST3 ) 25,92,25 
25 IF (CASE) 22,23,22 

22 Cl ( 78 )=C2 (78 ) + ( Cl (77 )— 077 ) * ( C278-CB278 ) / ( C177-CB17 7 ) 

C ( 78 )=C1 (78 ) 

C1(79)=C2(79)+(C1 (77)-C177)*(C279-CB279)/(C177-CB177> 

C ( 79 ) =C1 ( 79 ) 

GO TO 11 

92 CONTINUE 
WRITE (6,103) 

103 FORMAT ( 16HOCONVERGED STATE/25H0**************#"***#****#71H1 ) 
SAVPAP=1» 

IPRNT=NPR 
T£ST4=1 • 

ADINC=C2(79)-ALD1 
TATINC=TATD/Z-TAT1 
ALI NC=C2 ( 78 ) — ALFl 
ALT INC=ALTD-ALT1 

23 CASE=1. 

C ( 79 ) =C2 ( 79 ) 

C { 78 > =C2 ( 78 ) 

Cl (79)=C2 (79) 

Cl( 78 ) =C2 ( 78 ) 

GO TO 11 
END 


33 


SIBFTC FCALC 

CFCALC SUBROUTINE TO RUN TRAJECTORIES 
SUBROUTINE FCALC 

DIMENSION T( 150), C( 150) , ALTA 1250) >RHOTA( 250) »PRETA(250)»FVSTA(250) 
1 * FMT A ( 30 ) »CDTA { 30 ) *CLTA(30) 

DIMENSION Till 50) » Cl ( 150 ) 

DIMENSION C21150) ,T2( 150) 

COMMON T ♦ C 

COMMON ALTA»RHOTA*PRETA*FVSTA,FMTA»CDTA*CLTA 
COMMON C1,T1,TEST1,TEST2,C2 ,T2 

EQUIVALENCE (T1ME,T(2) ) , (STEP » T ( 3 ) ),(X»T(4)),(R,T(5))»(V»T(o))» 

1 ( TATER *T (7 ) ) * (SAM,T(8 ) ) * (ALAM2*T (9 ) ) , ( ALAM3 , T ( 10 ) > > ( ALAM4 » T ( 1 1 ) ) , 

2 ( ALAM5.K12 ) ) , (XDOT,T(13) ) , (ROOT. T( 14) ), ( VDOT»T( 15) ) * (TATDOT » 

3 T ( 1 6 ) ) > ( SAMDOT »T (17) ) , (DLAM2,T(18 ) ) , (DLAM3,T(19) ) , ( DLAM4 , T ( 20 ) J , 

4 (DLAM5.TI21) ) 

EQUIVALENCE (N»C(1 ) ) , (STEPP ,C( 2 ) ), (KA1P,C(3) ) , ( STEPC * C ( 4 ) ) > (K.A1C»C 
1(5) ) » ( A2 » C ( 6 ) ) * (A3>C(7) ) *(A4>C(8) ) » ( A5 > C ( 9 ) )* ( A6» C 1 10 ) )» IA7»C(11) ) 
2 

3» (60, Cl 13) ) » (GRAVO,C( 14) ) , (RO,C( 15 ) ) , (PRES0»C(16) ),( OMEGA ,C ( 17 ) ) ,( 
4PH I PR »C ( 1 8 ) ) » (AZ»C( 19) ) 

5 

6 > ( SAVPAP * C ( 2 1 ) ) » ( IPRNT »C( 22 ) ) , (NPR»C( 23) ) , ( I OPT ,C ( 24 ) ) ,(STP,C(25)) 

7 

8 j ( T IMPC »C ( 28 ) ) *CW0,C(29) ) » (WDO»C(30) ) »(FO,C(31) >, ( AEO » C ( 32 ) ) , ( AREA 

9 » C ( 33 ) ) j ( T 1M0 » C ( 34 ) ) 

EQUIVALENCE ( T IM1 ,C( 35 > ) » ( T IM2 ,C (36 ) ) , ( T IM3 ,C ( 37 >) » ( TI M4»C ( 38 ) > * ( T 
1IMF1 »C ( 39 ) ) , ( T IMF2 »C ( 40 ) ) » ( TIMF3 »C( 41 ) ) » ( T IMF4.CI 42 > ) > (FI ,C(43 ) ) , ( 
2F2,C(44))*(F3»C(45) )*(F4>C(46) )»(AE1»C(47) ),(AE2.C(48) ).(AE3»C(49) 

3 ) > ( AE4 » C ( 50 ) ), (TIMW1.C151) ) , ( TIMW2,C< 52 ) ) , (T1MW3,C(53) ),(TIMW4»C(5 
44) ),(TIMW5»C(55)),(TIMW6»C(56) ) , ( T I MW7 , C < 57 ) ) , ( WD1 »C ( 58 ) ) , (WD2»C(5 
59) ) , (WD3,C(60) ) >( WD4,C(61 ) ) , (WD5»C(62) ) , (WD6,C(e>3) ) > (WD7,C(64) ) 

6 

7 » ( ICUT » C ( 66 ) ) » ( TIMCC» C ( 67 ) ) , (TCL*C(68) ) » < Z1 »C (69) ) » (22 *C( 70) ) » ( Z3 » 
8C( 71) ) > (Z4 »C ( 72 ) ) > ( Z5 > C ( 73 ) ) , (Z6»C(74) ) ♦ <R2*C(75) ) 

9 

EQUIVALENCE (ALPHAP,C( 76) ) » ( ALFMAX ,C ( 77 ) ) , ( ALPHAC* C( 78 ) ) » ( DALF »C ( 7 
19) >, <DLAMX,C(80) ) » ( DELALF»C(81 ) ) » ( DELADT»C (82 ) ) 

2 

3 » ( WC »C ( 84 ) )»(SIP»C(85) ) »(FC»C(86)) 

4 

EQUIVALENCE ( ALPHA ,C ( 90 ) ) » ( SA, C ( 91 ) ) , { CA,C (92 ) ) , ( ST , C ( 93 ) ) * ( CT »C ( 9 
14) ) * ( SALT *C( 95 ) ) , ( CALT , C i 96 ) ), (GRAV,C(97) ) * ( ALT »C ( 98 ) ) »(DRAG,C(99) 
2 ) * ( FVS , C( 100) ) t (FMACH.Ct 101 ) ) » ( VS »C ( 102 ) ) * ( CD,C ( 103 ) ) » (QD *0(104) ) * 
3 (RHO*C( 105 ) ) * (PRES *C( 106) ) , ( AE »C ( 107) > » (PD IFFtC ( 108 ) ) » ( FRAC »C < 109 ) 

4) » (WGT.CC110) > * (WDOT ,C ( 111 ) ) » ( F ,C ( 1 12 ) ) » ( A ,C( 113) ) » ( PH I ,C ( 114 ) ) ♦ (C 
5HI »C( 115 > ) » ( ISAVE1 » C( 116) ) * ( ISAVE2 »C( 117) ) » ( JJ » C( 118 ) ) *( IMODE*C( 11 
69) >»(U,C( 120) > , (KA1»C(121 ) ) ,(CL,C(122 )) , (ALIFT»C( 123) ) 

7 

8 » ( TATD* C( 130) ) * (ALTD»C(131) > » ( T0L1 , C ( 132 ) ) , < T0L2, C ( 133 ) ) , ( T0L3 »C( 1 
934) ) 

IF (TEST1) 4,5,4 
5 IF (TEST2) 4,8,4 
8 CONTINUE 
ITEST2=1 
STEP-STEPP 
KA1=KA1P 
WGT =W0 
SAM=WGT/G0 
TIME=TIMO 
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WDOT =— WDO 

AE=AEO 

ALAM2=0.0 

AL AM3=0 .0 

ALAM4=0.0 

ALAM5=0.0 

DLAM2=0.0 

DLAM3=0.0 

DLAM4=0.0 

DLAM5=0.0 

IF (SAVPAP) 7,7,6 

6 WRITE ( 6 » 100 ) 

100 FORMAT (47H06IVEN AND PRECOMPUTED VALUES FOR LIFTOFF PHASE) 

7 IM0DE=1 
CALL PRESET 
DO 13 1=1,150 
C2i D=cm 

13 T2< I )=T(I I 

4 CONTINUE 

DO 14 1=1,150 
C( I )=C2 ( I ) 

14 T ( I ) =T2 ( I ) 

IM0DE=2 

VELS=SQRT(V*V+2.0*OMEGA*RO*V*COS <PHIPR)*SIN(TATER)-*SIN(AZ )+ OM 

1EGA*OMEGA*RO*RO*COS ( PH I PR )*COS ( PHI PR ) ) 

THETAS=ATAN(SQRT( 1.0-V*V*COSCTATER)*COS(TATER )/ (VELS*VELS ) )/ (V 

1#C0S ( TATER ) /VELS ) ) 

V=VELS 

TATER=THETAS 

ALPHA=ALPHAC 

STEP=STEPC 

KA1=KA1C 

F=FC 

WGT=WC 

SAM=WGT/GO 

TIMO=T IME 

SA=SIN(ALPHA) 

CA=COS (ALPHA ) 

ST=S IN ( TATER ) 

CT=COS ( TATER ) 

SALT=S1N< ALPHA+TATER) 

CALT=C0S ( ALPHA+TATER ) 

FRAC=RO/R 

H=V*V 

WDOT =— F/S IP 

SAMDOT=WDOT/GO 

VDOT =F/SAM # CA -GRAV *CT 

TATDOT=F/ ( SAM*V) *SA+ST*GRAV/V-ST*V/R 

B= (H*SALT )/ (V*DALF+6RAV*ST+SA*F/SAM-SA*H/R*CALT ) 

IF (B) 11,12,12 

11 ALAM2=-1 • 0 
GO TO 15 

12 ALAM 2 = 1 • 0 

15 ALAM3=ALAM2*B*CA/V 
ALAM4=ALAM2*B*SA 
ALAM5=0.0 

DLAM2=-2.0*ALAM3*GRAV/R*CT+ALAM4/R*ST*2.0*GRAV/V-ALAM4/R*ST*V/R 
DLAM3=-ALAM2*CT+ALAM4/V* ( t F*SA ) / ( SAM*V ) +1 GRAV/V+V/R ) *ST ) 
DLAM4=ALAM2*V*ST-ALAI13*GRAV*ST-ALAM4*GRAV/V*CT+ALAM4*V/R*CT 
DLAM5=F/SAM**2*< ALAM3*CA+ALAM4/V*SA ) 
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IF (SAVPAP) 1 0 > 10 1 9 
9 WRITE (6,200 

2o0 FORMAT (62H0G1VEN AND PRECOMPUTED VALUES FOR CALCULUS OF VARIATION 
IS PHASE) 

10 CALL CALVAR 
WGT=SAM*GO 
ALT-R-RO 
RETURN 
END 
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SIBFTC PRESE 

CPRESET SUBROUTINE FOR PRE SET ALPHA SEGMENT 
SUBROUTINE PRESET 

DIMENSION T ( 150) »C ( 150 ), ALTAI 250 ) >RHOTA<250) *PRETA(250> »FVSTA(250) 
1 » FMTA< 30 ) *CDTA( 30) »CLTA ( 30) 

DIMENSION C1(150)*T1( 150) 

DIMENSION C2(150),T2( 150) 

COMMON T » C 

COMMON ALTA,RHOTA»PRETA>FVSTA,FMTA,CDTA,CLTA 
COMMON Cl * T1 » TEST l>TtST2>C2»T2 

EQUIVALENCE I T I ME , T ( 2 ) ) » tSTEP,T(3) ) ,<X»T(4)),(R*T(5))*(V*T(6)), 

1 (TATERtT (7) )» ( SAM . T ( 8 >) * ( ALAM2 » T < 9 ) I » (ALAM3»T( 10) ) * l ALAM4 » T ( 11 ) ) » 

2 ( ALAM5,T (12 ) ) > <XD0T»TC13I ) » (RDOT,T ( 14) ) , ( VDOT,T( 15) ) » (TAT DOT » 

3 T ( 16 ) ) » ( SAM DOT » T ( 1 7 J ) , (0LAM2.T ( 18 > ) » ( DLAM3 » T ( 1 9 ) ) > ( DLAM4 > T ( 2 0 ) ) , 

A (DLAM5 »T (21 S ) 

EQUIVALENCE ( N , C (1 ) ) , ( STEPP * C ( 2 ) ) , (KA1P,C(3) ) , (STEPC»C( A) ) . (KA1C.C 
1(5) ) * ( A2»C(6) ) » (A3 »C( 7) ) » ( A4>C (8) ) , (A5 »C(9) ) » CA6»C (10) ) ♦ ( A7»C( 11 ) ) 
2 

3* (GO »C(13 ) ) » (GRAVO >C( 14 ) ) * ( RO »C ( 15 ) ) » (PRES0»C(16) ) ♦ (0MEGA»C<17) ) * ( 
APHIPR.CU8) ) > ( AZ»C(19) ) 

5 

6 » t SAVPAP * C < 2 1 ) > » (IPRNT,C(22) ) » <NPR,C<23) ) >(I0PT,C(24) ) , (STP»C(25) > 
7 

8* (TIMPC»C(28 ) 3 » (WO * C < 29 ) ) » (WD0»C(30 > ) » ( FO »C (31 > ) ♦ (AE0»C(32) ) , (AREA 
9>C«33) ) » ( T I MO »C (34 ) ) 

EQUIVALENCE ( T I Ml , C ( 3 5 ) ) * ( T IM2 ,C ( 36 ) ) » ( T I M3 ,C C 37 ) ) ♦ ( T IMA, C ( 38 ) ) , ( T 
1IMF1 »C(39 ) ) , (TIMF2*C(40) ) »(TIMF3»C( 41) )»(TIMEA»C(A2) )*(F1»C(43) )»( 
2F2»C(4A) ) , (F3.CU5 ) > , (F4»C(46) ) » (AE1,C(47) ) » ( AE2*C(48) ) . (AE3,C(49> 

3) » ( AE4,C( 50) ) • ( TIMWl.Ct 51 ) )*(TIMa2>C<52)) . (TIMW3.C (53 ) ) * (TIMW4»C(5 
44 ) ) > ( T I MW 5 » C ( 5 5 ) ) » (TIMW6>C(56) ) , (TIMW7»C(57) ) » ( WD1 »C ( 58 ) ) » (WD2*C(5 
59) ) * (WD3,C(60) ) » (WD4,C(61) ) » (WD5»C( 62) ) , (WD6»C(63) ) , <WD7»C(64) ) 

6 

7 > ( I CUT > C ( 66 ) ) > ( TIMCC,C(67) ) »( TOL»C(68) ) * ( Z1 *C ( 69 ) ) » ( Z2 * C ( 70 > ) » C 23 , 
8C(71) ) , (Z4,C(72) ), (Z5,C(73) ), (Z6»C( 74) ) , (R2*C(75) ) 

9 

EQUIVALENCE ( ALPHAP »C ( 76 ) ) , { ALFMAX , C ( 77 ) ) , (ALPHAC»C ( 78 ) ) , (DALF»C(7 
19) ) »<DLAMX,C(80>>, < DELALF, C 1 8 1 ) ) , ( DELADT ,C ( 82 ) ) 

2 

3»(WCsC(84))»(SIP»C(85))»(FC»C(86)) 

4 

EQUIVALENCE ( ALPHA »C ( 90 ) ) . ( SA , C ( 91 ) ) , ( CA.C ( 92 ) ) , ( ST , C ( 93 )) . ( C T » C ( 9 
14) ) » < SALT »C ( 95 ) ) > ( CALT, C ( 96 ) ) , (GRAV»C(97) ) , (ALT »C(98) ) , (DRAG, CC 99) 
2) » ( FVS > C ( 100) ) , < FMACM »C ( 101 ) ) » < VS »C ( 102 ) ) » (CD»C( 103) ) , ( QD.Ct 104 ) ) » 
3 ( RHO » C ( 105 ) ) ♦ (PRES, C( 106) ) » ( AE ,C ( 107 ) ) » ( PDI FF* C ( 108 ) ) » ( FRAC,C 1 109 ) 

4 ) ♦ ( WGT ,C( 110) > * (WDOT»C( 111 ) )j(r,C(112)),<A»C(ll3))»(PHI,C(114>)»(C 
5HIfC(ll5) )»( ISAVE1 1 C( 116) ) »< ISAVE2 >C( 117) ) » ( JJ*C( 118 ) ) » (IMODEiCdl 
69) ) » (U,C(120) > » (KA1,C(121) ) > ( CL. C 1122) ) » ( ALIFT.C( 123 ) ) 

7 

8 ? ( T ATD j C ( 130) ), ( ALTO, C( 131 )),( TOU * C( 132 ) ) » ( T0L2 , C ( 133 ) )»(T0L3»C(1 
934) ) 

ISAVE1=1 

ISAVE2=1 

CH0P=0. 

JJ=1 
IW=0 
I A= 0 

1ST AGE = 0 
TEMP=STEP 
KK= 1 



IPRINT=IPRNT-1 
IF ( T I Ml ) 60 >60 >6 1 

60 TMAX=T IMPC 
GO TO 62 

61 TMAX=T IM1 

62 IF (TIMF1) 63,63.64 

63 TMAXF=TIMPC 
60 TO 65 

64 TMAXF=T IMF1 

65 IF { TIMW1 ) 66.66,67 

66 TMAXW=T I MPC 
GO TO 6 

67 TMAXW=TIMW1 
GO TO 6 

7 GO TO <10.11.12) »KK 

10 JJ=JJ+1 

IF ( JJ-4 ) 51,51.22 

51 IF ( C ( JJ+34 ) ) 22.22,21 

22 TMAX=TIMPC 
GO TO 6 

21 TMAX=C( JJ+34) 

GO TO 6 

11 I A=IA+1 
AE-C ( I A+46 ) 

FO=C ( I A+42 ) 

IF ( IA-3) 52.52.24 

52 IF <C< IA+39) ) 24,24.23 

24 TMAXF=TIMPC 
I A=3 

GO TO 6 

23 TMAXF=C< IA+39) 

GO TO 6 

12 IW= IW+1 

WDOT =-C < IW+57) 

IF f IW-6) 53,53,26 

53 IF ( C ( I W+51 ) ) 26,26,25 
26 TMAXW=T 1 MPC 

I W=6 

GO TO 6 

25 TMAXW=Ct IW+51) 

6 IF ( TMAX.GT .TMAXF. OR.TMAX .GT.TMAXW ) IF ( TMAXF-TMAXW) 8,9.9 
TCHEK=TMAX 
KK = 1 

GO TO 13 

8 TCHEK=TMAXF 
KK=2 

GO TO 13 

9 TCHEK=TMAXW 
KK.= 3 

13 STEP=TEMP 
CH0P=0. 

1ST AGE= 1 

1 CALL INTG <T,N»K.A1,A2»A3,A4,A5,A6»A7) 

2 IF tCHOP) 40,32,33 
40 TIME=TCHEK 

GO TO 7 
33 CHOP=— 1 • 

GO TO 14 
32 CONTINUE 

TOGO=TCHEK“TIME 
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IF (TOGO-STEP) 30,14,14 
30 IF (TOGO) 5,7,5 
5 CONTINUE 
STEP=TOGO 
CH0P=1. 

GO TO 1 

14 IF (SAVPAP) 4,4,15 

15 IPRINT=1PRINT+1 

IF ( IPRINT-IPRNT ) 17,19,17 

17 IF (ISTAGE) 4»4»16 

19 IPRINT=0 

16 CALL OWT 
ISTAG£=0 

4 CALL INTM 

IF (AbS(TIME-TIMPC)-.OOOl) 18,2,2 

18 IF (SAVPAP) 20,20,29 
29 CALL OWT 

20 RETURN 
END 
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SIBFTC CALVA 

CCALVAR SUBROUTINE TO RUN CALCULUS OF VARIATIONS TRAJECTORY 
SUBROUTINE CAL VAR 

DIMENSION T C 15 0 ) »C ( 150 ) , ALTA ( 2 50 ) .RHOTA ( 250 ) ,PRETA(250 ) .FVSTAl 250 ) 
1 * FMT A I 30 ) *CDTA( 30) >CLTA( 30) 

DIMENSION Cl ( 1 50 ) »T1 l 150 ) 

DIMENSION C2(15O)»T2(150) 

COMMON T > C 

COMMON ALTA»RHOTA*PRETA,FVSTA»FMTA,CDTA,CLTA 
COMMON C1»T1,TEST1,T£ST2>C2»T2 

EQUIVALENCE <TIME»T<2 ) ) ,(STEP,T<3) ) ,(X>TC4)),(R>T<5))>(V>T(6)), 

1 l TATER >T (7 ) ! » (SAM>T(8) ) » ( ALAM2 *T ( 9 ) ) * ( ALAM3> T ( 10 J ) , (ALAM4,T(11 > ) , 

2 t ALAM5 » T ( 12 ) ) >(XD0T»T(13) I , (RD0T»T(14> ) » I VDOT , T 1 15 ) ) » ( TATDOT » 

3 T( 16! ) » (SAMDOT.Tt 17) ) , (DLAM2,T(18) ) . (DLAM3,T< 19) ) * ( DLAM4. T < 20 ) ) » 

4 (DLAM5 *T( 21 ) ) 

EQUIVALENCE (N»C(1) ) , ( STEPP * C ( 2 ) ),(KA1P»C<3) > , ( STEPC. C ( 4) ) »<KA1C*C 
1(5) )* (A2>C(6) ) » l A3 > C ( 7 ) ) > (A4>C(8) ) * ( A5 > C ( 9 ) ) > (A6»CtlO) > * ( A7 » C( 11) ) 
2 

3 > ( GO * C ( 13 ) ) > (GRAVO»Ct 14) > » (R0.CC15) ) , (PRES0.CU6) ) » ( OMEGA, C( 17) ) » ( 
4PHIPR*C(18) ) *( AZ»C(19) ) 

5 

6 » ( SAVPAP» C ( 21 ) ) > (IPRNT»C(22) ) , (NPR*C( 23 ) ) > ( IOPT »C ( 24) ) »(STP»C(25) ) 
7 

8 * ( TIMPC iC ( 28 ) ) » (W0*C(29) ) » (WD0»C(30) ) » ( FO * C ( 3 1 ) 1 » { AE0»C(32) ) ♦ ( AREA 
9»C(33) ) »(TIM0*C(34) ) 

EQUIVALENCE ( T IM1 , C ( 35 ) ) » ( T IM2 ,C ( 36 ) ) , { T IM3 , C ( 37 ) ) * < T IM4 »C ( 38 ) ) » ( T 
1 IMF1 >C ( 39 ) 3 » ( T 1MF2 »C( 40 ) )» (TIMF3»C(41) ) » { TIMF4»C< 42 ) ) » ( FI *C ( 43 ) ) » ( 
2F2,C(44) ) > (F3>C<45») » (F4»C(46) ) » (AE1»C(47) ) » ( AE2 »C (48 ) ) » ( AE3 »C( 49 ) 

3) * (AE4,C(50) )» ( T IMW1 > C ( 51 ) >, (TIMW2*C(52) ) » ( T1MW3»C ( 53 ) >»(T1MW4»C(5 
44) ) , ( TIMW5,C(55) ) , ( T I MW6 » C ( 56 ) ) , ( T IMW7 »C ( 57 ) ) » IWD1,C(58) ) » (WD2»C(5 
59) ) » (hD3>C(60) ) » <WD4»C(61) ) » (WD5.C( 62 ) ) ,(WD6»C(63) ) , (UD7»C(64) ) 

6 

7»(ICUT»C(66) )> ( T IMCC> C ( 67 ) ) , (T0L.C168) ) , (Z1,C ( 69) ) . (22 »C( 70) ) » ( Z3 » 
8CC71) ),(Z4.C(72) ),(Z5»C<73) ) >(Z6»C(74)) ,(R2*C(75I ) 

9 

EQUIVALENCE ( ALPHAP »C ( 76 1) » ( ALFMAX , C ( 77 ) ) , ( ALPHAC»C ( 78 ) ) » ( DALF »C ( 7 
19) ) , (DLAMX,C(80) ) , { DELALF,C( 81 ) ) , < DELADT » C < 82 ) ) 

2 

3»(WC»C(84) ) » (SIP»C(85) ) , (FC.C(86) 3 

4 

EQUIVALENCE ( ALPHA, C( 90) > » (SA»C{9l) ) , (CA»C(92) ) , ( ST >C i 93 ) ) ♦ ( CT >C ( 9 
14) ) » (SALT*C(95 > ) » (CALT*C(96) ) » (GRAV »C (97) ) »(ALT»C(98) ) » (DRA6»C(99 ) 
2) » (FVS»C(100) ) ♦ (FMACH*C(101) ) » ( VS»C ( 102 ) ) » (CD»C(103) ) » (QD»C( 104) ) , 
3 ( RHO> C ( 105 ) ) , (PRES*C( 106 ) ) » ( AE »C ( 107 ) ) » ( PD IFF, C ( 108 ) ) * ( FRAC »C ( 109 ) 

4) > < WGT >C( 110 ) ) » (WDOT >C( 111 ) ) » < F.C ( 1 12 ) ) > C A,C( 113) > > ( PHI ,C< 114 ) ) » < C 
5HI,C(115))»USAVEi,C(116) ) » ( ISAVE2>C( 117) ) > ( J J > C < 1 18 ) > *(IM0DE»C(11 
69) ) > ( U»C( 120 ) ) » ( KA1 »C ( 121 ) )»(CL»C(122))»(ALIFT»C(123) ) 

7 

8* ( TATD»C( 130 ) ) , < ALTD,C( 131) ) » ( TOLl , C ( 132 ) ) , ( T0L2 ( 133 ) ) , (T0L3»C(1 
934) ) 

I PRINT= IPRNT 
INISH=IPRINT— 1 
SEARCH=1*0 
CUTOFF=TIMCC 
JET=1 

TLAST=TIME~1. 

1 CONTINUE 

GO TO ( 31 »32 ) > I CUT 
32 CONTINUE 
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VALUE=TIME 
RATE= 1 • 0 
GO TO 33 
31 CONTINUE 
VALUED 
RATE^VDOT 

CUT0FF=Z1*SQRT( t (1 . 0+Z2 ) *Z3* ( 1 *0+ ( Z4*R ) /R2 ) )/( ( 1 *0+ ( Z4*R*R ) / { R2*R2 
1)*ST*ST)*R)+Z5)+Z6 

33 CONTINUE 
9 CONTINUE 

IF < ABS (VALUE -CUTOFF ) -TOL ) 11,12,1 2 

12 IF (SEARCH) 14,14,13 

13 GO TO (15,16) , JET 

15 CALL INTO < T ,N , K.A1 » A2 , A3 , A4 , A5 , A6 , A7 } 

WHI CH= ( VALUE-CUTOFF ) / ABS ( VALUE -CUTOFF ) 

JET=2 

16 CHANGE- ( VALUE“CUTOFF ) /ABS ( V ALU E^CU TOFF ) 

I F ( CHANGE / WH I CH ) 1 4 , 1 7 , 1 7 

14 STEP=(CUT OFF-VALUE ) /RATE 
SEARCH=0.0 

KA1 = 1 
GO TO 15 

17 IF (SAVPAP) 18,18,42 

42 IF ( ( T I ME- TL AST ) ^SEARCH ) 18,18,43 

43 INI SH= INI SH+1 

IF ( IPRINT-INISH) 18,19,18 
19 CALL OWT 
INI SH=0 

18 TLAST=T IME 
CALL INTM 
GO TO 1 

11 IF (SAVPAP) 35,35,34 

34 CALL OWT 

35 CONTINUE 
RETURN 
END 
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SIBFTC DAUX 

CDAUX SUBROUTINE TO EVALUATE DERIVATIVES 
SUBROUTINE DAUX 

DIMENSION T< 150 ) »C( 150 ) , ALTA ( 2 50 ) »RH0TA ( 250 ) , PRET A ( 250 ) >FVSTA ( 250 ) 
1 » FMTA ( 30 ) tCDTAC 30) ,CLTA(30) 

DIMENSION Cl ( 150 1 * T 1 ( 150 ) 

DIMENSION C2 ( 150 ) > T2 l 150 ) 

COMMON T » C 

COMMON ALTA jRHOTA, PRETA.FVSTA, FMTA »CDTA,CLTA 
COMMON C1»T1»TEST1,TEST2»C2»T2 

EQUIVALENCE CT1ME,T<2 ) ) , (STEP»T(3) ) »(X»T(4)),<R,T<5>>»(V.T<o)), 

1 ( T ATER »T t 7 > > , ( SAM » T ( 8 ) > > < ALAM2,T<9) > , (ALAK3,T< 10) ) , ( ALAM4,T < 11 ) ) , 

2 < ALAM5 »T ( 12 ) ) » ( XDOT >T(13)),(RD0T»T(14) > , ( VDOT » T ( 15) ) * I TATDOT ♦ 

3 T ( 16 ) ) > I SAMDOT »T( 17 J ) , (DLAM2.TU8) ), <DLAM3,T(19) ) , <DLAM4,T<20 ) ) , 

4 ( DLAM5 »T ( 21 ) ) 

EQUIVALENCE [ N , C (1 ) ) , ( STEPP > C < 2 )) , ( KA1P >C ( 3 ) ) » ( STEPC , C ( 4 )) » ( KA1C ,C 
1(5) ) > (A2»C(6) ) , <A3,C<7) ) » ( A4,C(8) ) , ( A5,C(9)) » (A6»C( 10) }»(A7,C< 11) ) 
2 

3 » ( GO » C ( 13 ) ) * (6RAV0»C(14) ) » <R0»C(15) ) » < PRESO >C 1 1 6 ) )» (OMEGA, C( 17 > ) , ( 
4PHIPR ,C ( 1 8 ) ) ,(AZ,C(19) ) 

5 

6, ( S AVPAP » C ( 21 ) ) » (IPRNT,C(22) ), ( NPR »C ( 23 ) ) > ( I0P7 ,C<24> ) , (STP,C ( 25 > ) 
7 

8, ( T IMPC ,C ( 28 ) ) * (W0,C(29) ) » IWDO*C(30) ) , ( FO »C ( 3 1 ) ) » (AEO»C(32 ) ) , (AREA 
9 » C ( 33 ) ) * ( T I MO >C 1 34 ) ) 

EQUIVALENCE ( T IM1 ,C( 35 ) ) » ( TIM2 ,C ( 36) ) , ( T IM3 ,C <37 ) ) » ( T IM4 ,C ( 38 ) > » ( T 
II MF1 » C ( 39 ) ) , (TIMF2»C(40) ) » (TIMF3,C(41) > » ( T IMF4, C ( 42 ) ) , < FI ,C ( 43 ) ) , ( 
2F2,C(44) ) , (F3,C(45 ) ) , ( F4,C(46) ) , (AE1,C(47) ) , l AE2»CC48) ) » <AE3»C(49) 

3) , ( AE4,C( 50) ) , (TIMW1,C(51) ) * (TIMW2,C(52) ) , ( TI MW3, C ( 53 ) ),(TIMW4,C(5 
44) ) , (TIMW5,C(55) > , ( TIMW6 ,C( 56) ) , <T IMW7 »C< 57 ) ) , ( WD1.CC58) ) , (WD2,C(5 
59))»(WD3»C160)),( WD4»C ( 61 ) ) » ( WD5 >C ( 62 ) ) » ( WD6 , C ( 63 ) ) , (WD7,C(64) ) 

6 

7* ( I CUT »C ( 66 ) ) . (TIMCC,C(67) ) , (T0L,C(68 ) ) , ( Z1 ,C ( 69 ) ) » ( Z2 ,C ( 70 ) ) » ( 23 » 
8C(7D ) » (Z4,C(72) ) , (Z5»C(73) ) ,<Z6,C<74) ) , <R2,C<75) ) 

9 

EQUIVALENCE (ALPHAP ,C( 76 ) ) , ( ALFMAX ,C ( 77 ) ) , ( ALPHAC * C ( 78 ) ) , (DALF,C(7 
19 ) ) , (0LAMX»C(80) ) , ( DELALF » C ( 81 ) ) , ( DELAUT,C(82 ) ) 

2 

3, ( WC » C ( 84 ) ) , (SIP,C(85) ) , (FC»C(86) ) 

4 

EQUIVALENCE { ALPHA, C( 90) ),(SA,C(91) ),(CA,C(92)) » ( ST»C (93 ) ) » ( CT ,C ( 9 
14) ) , (SALT ,0 95) ) »<CALT,C<96) ) , (GRAV,C<97) ) , <ALT»C<98) ) , (DRAG,C(99) 
2) » ( FVS.Ct 100) ) , (FMACH,C( 101 ) ) , ( V$»C ( 102 > ) » ( CD »C( 103 ) ) , <GD,C( 104 ) ) , 
3(RH0,C( 105 ) ) , (PRES,C( 106 ) ) » ( AE »C ( 107 ) ) , < PDI FF , C ( 108 ) ) , (rRAC,C( 109) 

4) , ( WGT,C(11C) ) , (WDOT»C< 111) ),(F,C(112)),(A»C(113))»(PHI ,C{ 114) ) , (C 
5H I , C ( 115 ) ) » ( ISAVE1 »C ( 1 16 ) ) » ( ISAVE2 »C( 117 ) ) , ( JJ»C( 118) ) , ( I MODE, C ( 11 
69) ) » (U,C( 120) ) , (KA1,C< 121) ) , (CL,C( 122 ) ) ,(ALIFT,C( 123) > 

7 

8, (TATD»C( 130) ) » ( ALTD ,C { 1 31 ) ) , ( T0L1 ,C( 132 ) ) * (T0L2»C(133 ) ) » ( T0L3,C<1 
934) ) 

GO TO (1,2), IMODE 
1 CONTINUE 

GO TO (13,20,14,30,13) »JJ 

13 ALPHA=ALPHAP 
GO TO 3 

20 A LPHA=ALPHAP+( ALFMAX— ALPHAP ) *( T I ME- T I Ml )/(TIM2-TIMl) 

GO TO 3 

14 ALPHA=ALFMAX 
GO TO 3 
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30 ALPHA=ALPHAP+( ALFMAX-ALPHAP )#( TIM4-TIME ) / ( TIM4-T1M3) 

3 SA=S IN (ALPHA) 

CA=C0S ( ALPHA) 

GO TO A 
2 CONTINUE 

U=SQRT ( ALAM4**2+V*V*ALAM3**2 ) 

SA=ALAM4/U 

CA=ALAM3*V/U 

ALPHA=ATAN2(SA,CA) 

A ST=S IN ( TATER ) 

CT=COS(TATER) 

FRAC=RO/R 

GRAV=GRAVO*FRAC*FRAC 
XDOT=FRAC*V*ST 
RDOT =V*CT 

GO TO (5,9) * IMODE 

5 CONTINUE 
ALT=R-RO 
SAMDOT=WDOT /GO 
CALL TABLE 

PD I FF=PRESO-PRES 
F=FO+AE*PDI FF 

VDOT= ( F *CA-DRAG ) /SAM-GRAV*CT 
IF (V) 7, 6, 7 

6 TATDOT=0»0 
60 TO 8 

7 TATDOT=(F*SA+ALIFT )/<SAM*V)+(GRAV/V-V/R)*ST 

8 RETURN 

9 SAMDOT =—F/S IP/GO 

VDOT =F/SAM * CA -GRAV *CT 
T ATDOT =F/ ( SAM*V )*SA+ST*GRAV/V-ST*V/R 

DLAM2=-2.0*ALAM3*GRAV/R*CT+ALAM4/R*ST*2.0*GRAV/V-ALAM4/R*ST*V/R 
DL AM3=-ALAM 2*CT+AL AM4/V* ( ( F*SA ) / ( SAM*V ) + ( GRAV /V+V/R ) *ST ) 
DLAM4=ALAM2*V*ST-ALAM3*GRAV*ST-ALAM4*GRAV/V*CT+ALAM4*V/R*CT 
DLAM5=F/SAM**2*(ALAM3*CA+ALAM4/V*SA) 

RETURN 

END 
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SIBFTC OWT 

COWT OUTPUT SUBROUTINE 
SUBROUTINE OWT 

DIMENSION T(150)»C(150) .ALTAI 250 ) ,RHOTA(250> »PRETA<250) »FVSTA(250) 
1 , FMTAI 30) »CDTA(30) .CLTA ( 30) 

DIMENSION 0(150] ,T1( 150) 

DIMENSION C2 ( 150 ) . T2 ( 150 ) 

COMMON T . C 

COMMON ALTA.RHOTA.PRLTA.FVSTA.FMTA.CDTA.CLTA 
COMMON C1.T1.TEST1.TEST2.C2.T2 

EQUIVALENCE (TIME. T < 2 > ) , ( STEP. T(3)),(X»T(4)).CR»T<5))»(V»T(6))» 

1 ( TATER »T ( 7 J ) . ( SAM.T (8 ) ) . (ALAM2 .T( 9 ) ) . ( ALAM3 » T ( 10 ) ) , < ALAM4,T( 11 ) > . 

2 ( ALAM5.T112)) . (XDOT.T( 13) ) , ( ROOT » T ( 14 ) ) , (VDOT.T115) ) , (TATDOT. 

3 T( 16)) » ( SAMDOT »T ( 17 ) ) . (DLAM2.T (18 ) ) . (DLAM3.T (19) > > ( DLAM4.T1 20) ) » 

A ( DLAM5.T ( 21 ) ) 

EQUIVALENCE (N,C(1 ) ) , CSTEPP,C(2) ) . <KA1P.C<3> ) . (STEPC.C(A) ) , (KA1C.C 
1(5) ) > (A2»C(p) ) > (A3»C(7)) . (AA.C18) )»(A5,C(9)) , CA6» C(10) )»(A7»C(11J ) 
2 

3 > ( GO »C ( 13 ) ) > (GRAV0.CU4) ) .<R0.C(15) ) > ( PRESO . C ( 16 ) ) . ( OMEGA, C( 17 ) ) . ( 
4PH1 PR , C ( 18 ) ) ,(AZ,C<19>) 

5 

6, ( SAVPAP . C ( 2 1 ) ) , ( IPRNT.C122) ) , (NPR»C(23) ) , (10PT »C(24) ) , (STP.C(25) ) 
7 

8. (TIMPC.C(28) > > (W0.C(29) ) . (WD0.C13O) ) ♦( F0.C(31) ).(AE0.C(32) ) . (AREA 
9»C(33) ) »(TIM0»C(34) ) 

EQUIVALENCE l T IM1 , C ( 35 ) ) . ( T IM2 , C ( 36 ) ) . ( T IM3 , C (37 ) ) . ( T IM4 . C ( 36 ) ) . ( T 
1IMF1 »C( 39 ) 3 > ( T IMF2 »C ( 40 ) ) , ( T IMF3 .C ( 41 ) ) » ( T IMF4. C (42 > ) , ( FI , C ( 43 ) ) , ( 
2F2 »C ( 44 ) ) , (F3»C(45) ) . (F4.C(46) ) , (AE1 ,C(47) ) , ( AE2 »C ( 48 ) ) , ( AE3 ♦ C ( 49 ) 

3) , (AE4.CI 50) ) , (TIMWl.Ct 51) ) » ( T I MW2 . C ( 52 ) > , (TIMW3.C(53) ) , (TIMW4.C15 
44) ) » ( TIMW5.C(55 ) ) » ( T 1MW6 »C ( 56 ) ) , ( T IMW7,C(57) ) . ( WDl >C(56) > . (WD2»C(5 
59 ) ) . (WD3»C(oO) ) » ( W04.C( 61 ) ) . ( WD5»C(62) ) ♦ CWD6»C(63) ) . (WD7»C(64) ) 

6 

7, ( I CUT, C( 66) ) , ( T1MCC.C ( 67 ) ) , (TOL,C(68) ),CZ1,C(69)),(Z2,C(70)),IZ3, 
8C (71) )♦ (Z4»C(72) ) ,(25.C(73) ) »(Z6»C(74) ) » (R2.C(75) ) 

9 

EQUIVALENCE ( ALPHAP , C ( 76 ) ) , ( ALFMAX »C(77 ) ) , (ALPHAC,C( 78) ) , (DALF»C(7 
19) ) » ( DLAMX »C (80) ) ,(DELALF»C(81) ) , ( DELADT »C ( 82 ) ) 

2 

3, (WC,C<84) ), (SIP,C(85J) , (FC,C(86) ) 

4 

EQUIVALENCE ( ALPHA ,C ( 90 ) ) » ( SA, C ( 91 ) ) » ( CA,C ( 92 ) > » ( ST »C ( 93 ) ) , ( CT , C < 9 
14) ) , l SALT »C ( 95 ) ) , ( CALT , C ( 96 ) ) » ( GRAV.C (97 ) ) , (ALT,C(98) ) , (DRAG,C(99) 
2) *<FVS»C< 100) ) , (FMACH.C(lOl) ) , ( VS ,C ( 102 ) ) , ( CD.C < 103 ) ) , (QD,C( 104 > ) , 
3( RHO,C( 105 ) ) » (PRES »C( 106 ) ) . ( AE , C C 107 ) ) . ( PD I FF, C ( 108 ) ) , ( FRAC »C ( 109 ) 

4) » (WGT.C(llO) ) . (WDOT.C(lll) ) ♦ ( F,C( 112 ) > , ( A.C( 113) ) » ( PHI ,C< 114 ) ) , (C 
5HI , C ( 115 ) ) » ( ISAVE1 ,C( 116 ) ) , ( ISAVL2 ,C( 117) ) , ( JJ, C( 1 18 ) ) , ( ImODE.CUI 
69) ) » (U,C( 120) J » ( KA1 , C ( 121 ) ) , ( CL ,C( 122 ) ) . ( AL I FT , C (123 ) ) 

7 

8, (TATD,C( 130) > » (ALTD.Ct 131 ) ) , ( T0L1 , C ( 132 ) ) , (T0L2,C( 133) ) , (T0L3,C(1 
934) ) 

Z=57. 295775 
E2=T ATER*Z 
E3 = T ATDOT*Z 
E4=ALPHA*Z 
PHI =X*Z/R0 
CH1=E2+E4+PHI 
WGT=$AM*GO 
GO TO (1,2) .IMODE 
1 CONTINUE 
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GO TO 13.4.3,5*3) , J J 

3 ADOT=0. 

GO TO 6 

4 AOOT= ( ALFMAX-ALPHAP ) / < T I M2-T IM1 ) 

GO TO 6 

5 ADOT=( ALFMAX-ALPHAP ) / < T I M3-TIM4 > 

GO TO 6 

2 ALT=R-RO 

DALF=ALAM2*V*SALT/U-GRAV*ST/V+SA*V*CALT/R-SA*F/ (V*SAM) 

ADOT=DALF 

A=ALAM2*RDOT+ALAM3*VDOT+ALAM4*TATDOT+ALAM5*SAMDOT 

6 ADOT=ADOT*Z' 

WRITE (6.100) X»XD0T»TIME»E4,AD0T,R»RD0T»ALT.WGT,WD0T,V,VD0T.F,SAM 
1 .SAMDOT » E2 > E3 »GRAV »PHI »CH I 

GO TO (11,12) .IMODE 

11 WRITE (6.200) RHO,PRES,VS,0RAG,ALIFT,O0,AE,FMACH.CD,CL 
GO TO 13 

12 WRITE (6.3C0) SI P, ALAM2, ALAM3 .ALAM4.ALAM5 » A.DLAM2.DLAM3 >DLAM4» 

1 DLAM5 

13 RETURN 

100 FORMAT (8H0XXX E14.7.12H XDOT E14.7.12H TIME E14.7, 

112H ALPHA £14.7. 12H ALDOT E14.7/8H RRR E14.7.12H 

2RD0T E14.7.12H ALT E14.7.12H wLBS E14.7.12H WOO 

3T E14.7/8H VEL E14.7.12H VDOT E14.7.12H THRST E14. 

47.12H MASS E14.7.12H MDOT E14.7/8H THETA E14.7.12H 

5 THDOT E14.7.12H GRAV E14.7.12H PHI E14.7.12H C 


6H1 E14.7) ^ 

200 FORMAT (8H RHO E14.7.12H 

H2h drag E14.7.12H 

2AEX E14.7.12H MACH 

3 E14.7 ) 

300 FORMAT (8H ISP L14.7.12H 

112H LAMB4 E14.7.12H 

2DLAM2 E14.7.12H 0LAM3 

3M5 E14.7) 

END 


PRES E14.7.12H 
LIFT E14.7/8H Q 
E14.7.12H CD 


VS E14.7, 

E14.7.12H 
E14.7 .12H CL 


LAMB2 E14.7.12H LAMB3 E14.7, 

LAMB5 E14.7/8H AAA E14.7.12H 
E14.7.12H DLAM4 E14.7.12H DLA 
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5IBFTC TABLE 

CT ABLE SUBROUTINE FOR LOOK UP AND AERODYNAMIC PARAMETERS 
SUBROUTINE TABLE 

DIMENSION T ( 150 ) ,C ( 150 )» ALTA ( 250) ,RHOTA ( 250 ) » PRETA< 250 ) ,FVSTA<250) 
1 > FMT A ( 30 ) »CDTA( 30 ) *CLTA ( 30 ) 

DIMENSION C1(150>,T1(150) 

DIMENSION C2(150),T2(l50) 

COMMON T > C 

COMMON ALTA* RHOT A*PRETA»FVSTA» FMT A , COT A » C LT A 
COMMON Cl » T1 » TEST1 »TEST2>C2*T2 

EQUIVALENCE ( T I ME,T ( 2 )) , ( STEP » T ( 3 ) ) , ( X»T ( 4 ) ) , (R , T ( 5 ) ) ♦ ( V»T ( 6 ) ) , 

1 (TATER,T(7) ) » (SAM»T (8 ) ) * ( ALAM2 * T < 9 ) ) * ( ALAM3* T ( 10 ) > » < ALAMA.T ( 11 > ) * 

2 1 ALAM5 * T ( 12 ) ) * ( XDOT , T < 13 } ) ,(RDOT,T( 14) ) * (VDOT.T (15) ) » ( TATDOT » 

3 T ( 16 ) ) , (SAMDOT»T(17) ) * (DLAM2»T(18 ) ) * { DLAM3 » T ( 1 9 ) ) * ( DLAM4 , T ( 20 > ) » 

4 ( DLAM5 * T ( 2 1 ) ) 

EQUIVALENCE ( N » C ( 1 ) ) , ( STEPP ,C( 2 ) ), ( K.A1P ,C ( 3 ) ) * ( STEPC *C ( 4 ) )»(KA1C,C 
l(5))»(A2»C<6))»(A3»C(7))»(A4»C(8))>(A5fC(9))*(A6*C(10))»(A7*C(ll)l 
2 

3 * ( GO »C { 13 ) ) » (GRAV0»C( 14) ) > (R0»C(15 ) ) . (PRESO,C( 16 ) ) * < OMEGA »C ( 17 ) ) , ( 
4PHIPR,C(18))*«AZ,C(19) ) 

5 

6 * (SAVPAP >C ( 21 ) ) * (IPRNT»C(22 ) ) > (NPR*C(23) ) » ( IOPT»C( 24) ) » <STP>C<25) ) 
7 

8* ( T IMPC * C C 28 ) ) * (W0,C(29 ) ) , (WD0tC(3O) ) » ( F0,C(31 ) ) * ( AE0.CI32 ) ) , ( AREA 
9, C( 33 ) )* ( T IMG * C ( 34 ) ) 

EQUIVALENCE ( T I Ml , C ( 35 ) ) ♦ ( T IM2 »C ( 3 6 ) ) • ( T I M3 ,C ( 37 ) ) » < T IM4*C t 38 ) ) * ( T 
1 1 MF1 * C ( 39 ) ) » { T I MF2 *C(40) > » ( TIMF3,C(41) ) , ( T IMF4 , C ( 42 ) ) * ( FI ,C (43 ) ) , ( 
2F2,C<44) ) , ( F3,C(45 > ) , ( F4,C< 46) ) , ( AE1 ,C<47 ) ) , < AE2,C(48 ) ) , ( AE3*C(49 ) 

3) * ! AE4,C(50) ) . ( TIMW1 *C < 51 ) ) * ( T IMW2 » C ( 52 ) ) , (TIMW3*C(53 ) ) » (TIMW4»C(5 
44) ) , (TIMW5,C<55 ) ) > ( T IMw6 » C ( 56 ) ) » ( T IMW7 *C( 57 ) ) , <WD1»C(58) > * (WD2*C(5 
59 ) ) * ( WD3 9 C ( 60 ) ) 9 ( VD4 9 C { 61 ) ) 9 (WD5 9C( 62) ) 9 ( WD6,C(63> ) 9 (WD7»C(64) ) 

6 

7 9 ( I CUT 9C( 06 ) ) 9 ( TIMCC*C(67) ) , (TOL*C(68) ) 9 ( Z1 9 C ( 69 ) ) 9 ( Z2 . C ( 70 ) ) 9 ( Z3 9 
8C(71 ) ) 9 (Z49C172) ) , (Z5*C(73) ) ,(Z6*C( 74) ) 9 (R2*C(75 ) ) 

9 

EQUIVALENCE ( ALPHAP , C ( 76 ) ) 9 (ALFMAX,C(77 ) ) , ( ALPHAC »C ( 7 8 ) ) , (DALF*C(7 
19) ) 9 ( DLAMX 9 C ( 80 ) > 9 ( DELALF 9 C ( 81 ) ) 9 ( DELADT 9 C ( 82 ) ) 

2 

39 (WC,C(84) ) 9<SIP»C(85) ),<FC,C(86>) 

4 

EQUIVALENCE { ALPHA 9 C ( 90 ) ) 9 ( SA 9 C { 91 ) > , ( CA, Ct 92 ) 1 9 ( ST 9 C ( 93 ) ) 9 ( CT 9 C ( 9 
14) ) 9 (SALT tC (95) ) ,(CALT,C(96 ) > , ( GRAV »C ( 97 ) ) ,(ALT*C(98) ) » ( DRAG ,C ( 99 ) 
2) 9 (FVS 9 CUOO) ) 9 (FMACH9C<101 ) > 9 (VS9CC102) ) ,(CD,C( 103) ) 9(QD9C(104 ) ) 9 
3(RH0,C(105) ), (PRES 9C( 106) 1 9 ( AE , C ( 107 ) ) , ( PDI FF ,C ( 108 ) ) » ( FRAC ♦ C ( 109 ) 

4) 9 (WGT.CdlO) ) 9 Ur DOT 9 C l 111) ),(F,C( 112) ) 9(A,C(113) )» (PH1,C(114) ) »(C 
5H1»C( 115) ) » ( ISAVtl * C ( 116) ) » ( ISAVE2 »Ct 117) ) » ( JJ *C( 118 ) ) * ( 1M0DE*C(11 
69) ) »(U.CU2C) ) *<KA1,C(121) ) ,(CL,C( 122) ) ♦ ( ALIFT *C( 123 ) ) 

7 

89 (TATD,C(130) ) * ( ALTD 9 C ( 131 ) ) . (T0L1»C(132) ) * ( T0L2 9 C ( 133 ) )»(T0L3 9 C (1 
934) ) 

GL0M=1. 0548599/GO 
I=ISAVE1 

3 IF (ALT) IOO 9 110, 110 
100 1=1 

IL = 1 
IU = 2 

GO TO 10 

110 IF ( ALT— ALTA ( I > ) It IO 9 4 

4 IF ( 1-250) 5,696 
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5 1=1+1 

go to no 

6 1=250 

I L=249 
IU=250 
GO TO 10 

7 1=1-1 

IF (ALT— ALTA( 1 3 ) 7,10,11 

10 FVS=FVSTA ( I ) 

RH0=RH0TA ( I )*GLOM 
PRES=PRETA(I) 

IF (1-250) 1,12,6 
1 IF (1-13 100,12,41 

41 I L= I — 1 
IU=I+1 
GO TO 12 

11 I L= I 
IU=I+1 

FRAC1= (ALT-ALTA (IL) >/( ALTAI IU) -ALTAI ID ) 
FVS s FRACl* ( FVSTAI IU)-FVSTA( ID )+FVSTA(lL) 

RHO= ( FRAC1* £ RHOT A ( IU)— RHOTA( IL ) }+RHOTA{ IL) )*GLOM 
PRES=FRAC1* ( PRETA ( IU ) —PRETA 1 1 L > ) +PRETA ( I L ) 

12 ISAV£l=l 
FMACH=V*FVS 
J=ISAVE2 

16 IF ( FMACH ) 200, 210, 210 
200 J=1 

JL= 1 
JU= 2 

GO TO 23 

210 IF ( FMACH-FMTA ( J ) ) 20, 23, 17 

17 IF (J-30) 18,19,19 

18 J=J+1 

GO TO 210 

19 J=30 
JL=29 
JU=30 

GO TO 23 

20 J=J-1 

IF (FMACH-FMTA(J) ) 20,23,24 

23 CD=CDT A ( J ) 

CL=CLTA( J ) 

IF (J-30) 14,25,19 
14 IF (J-l) 200,25,51 
51 JL=J— 1 
JU=J+1 
GO TO 25 

24 JL= J 
JU=J+1 
K=JL 

FRAC2= ( FMACH- FMT A(JL) )/(FMTA(JUJ -FMT A < JL) ) 
CD=FRAC2*(CDTA( JU)-CDTA( JL) )+CDTA( JL) 

CL=FRAC2* I CLT A( JU) -CLTA I JL ) ) +CLTA ( JL ) 

25 ISAVE2=J 
VS=1 •O/FVS 
QD= » 5*RH0*V*V 

DRAG=AREA*QD* ( CD*CA+ALPHA*CL*SA 5 
AL I FT= AR£A*QD* ( -CD*SA+ALPHA*CL*CA ) 

RETURN 

END 
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•5SIBFTC MATS 

SUBROUTINE MATS (A»X»N»M iNSING) 
DIMENSION A(5»6)»X(5.1) 

MM=N+M 
DO 15 1=2, N 
70 11=1-1 

7 DO 15 J = 1 , 1 1 

8 I F ( A ( I ,J) 19*15.9 

9 IF t ABS (A(J*J) )-ABS(A( I. J1 1)11.10.10 

10 R=A ( I » J 1 / A < J » J 1 
GO TO 130 

11 R=A I J .J 1 / A ( I.J) 

DO 12 K=1»MM 
B=A( J.K) 

A ( J » K 1 =A ( I ,K) 

12 At I »K)=B 
130 JJ=J+1 

13 DO 14 K=JJ.MM 

14 A(I»K)=A(I,K.)-R*A(J»K) 

15 CONTINUE 

IF ( ABStAtN.N) 1-1. E-081 16.16.17 

16 NS I NG=1 
RETURN 

17 DO 28 J=1 . M 
KK=N+ J 

X(N»J)=A(N»KK)/A(N.N) 

DO 28 1=2. N 
J J=N- 1+1 
B=0.0 
1 1 =N-I+2 
DO 25 K=I I .N 
25 B=B+A(JJ,K)*X(k.,J) 

I F ( ABStAt JJ. JJ1 1 — l.t-08 116.16.28 
28 Xt JJ.J)=(A(JJ.KK)-B)/AC JJ.UJ! 

NS I NG = 0 

RETURN 

END 
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SIBFTC ISO 
C1SO 

SUBROUTINE ISO t ul * G2 * Q3 , w4 ,05 * 06 * I 07 ) 

DIMENSION X(lS))»Y(15hu(5t6),C0O,l) # ICNT ( 3 } » R1 ( 3 } *R2 ( 3 } , DVAL ( 3 ) 

VAL=Q 1 

VARY-Q2 

VALD=Q3 

F I NC=Q4 

T0L=Q5 

T EST=Q6 

JG0=IQ7 

ICOUNT=ICNT ( JGO) 

K= 5* JGO 
K0=K-4 

IF ( ABS(VAL-VALD)-TOL) 34*54*4 

4 IF (TEST) 1,3*1 
3 I COUNT=0 

TEST=1.0 
M= 1 

R 1 ( JGO )=ABS(1*/VARY) 

R2UG0)=ABS< l./VAL) 

X ( KO ) =VAR Y^Rl ( JGO ) 

Y ( KO ) =VAL*R2 ( JGO ) 

VARY=X ( KO ) 

VAL=Y (KO) 

OVAL ( JGO ) =VALD*R2 ( JGO ) 

GO TO 53 

1 X(K)=X(K~1) 

Y(K)=Y(K-1 ) 

K = K-1 

IF (K-KO) 2,2*1 

2 X(K)=VARY*R1( JGO) 

Y ( K ) =VAL*R2 ( JGO ) 

K1=K+1 

I COUNT= I COUNT +1 
IF ( I COUNT-4 ) 6,6*5 

5 IC0UNT=4 

6 CONTINUE 
VARY=X(K) 

VAL= Y ( K. ) 

GO TO ( 12*23*25*31 ) *ICOUNT 
12 KANT=0 

11 VI NC= ( OVAL ( JGO) -VAL ) * ( VARY-X ( K1 ) ) / ( VAL-Y ( K1 ) ) 

VARY=VARY+VINC 
IF (KANT) 51,51,22 

22 IF ( ABS ( VI -VARY )-ABS < V2-VARY ) ) 16,16,17 

16 VARY=V1 
GO TO 52 

17 VAR Y =V2 

52 IF t ( VARY-X(K) )/VINC) 2 7*27*51 
27 VARY=X(K)+VI,NC 
GO TO 51 

23 N=3 

DO 24 I =1 * N 
L=K+ I “1 

Q( 1*1 )=X<L)**2 
Q ( I * 2 ) =X < L ) 

Q( I * 3 )«YCL) 

24 Q(I*4)-1*0 
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CALL MATS ( 0, CO, N, H , f’Sl N6 ) 
IF (NSING) 71.71.12 

71 CONTINUE 
A=CO ( 1.1) 

8 = 0 . 

c=o. 

D=CO 12*1) 

E=CO (3.1) 

KANT=1 
GO TO 19 

25 N =4 

DO 26 1=1, N 
L=K+I-1 

OU »1)=X(L)**2 
Q ( 1 »2>=X(L)*Y(L! 

Oil . 3 ) =Y ( L ) **2 
Q { I ,4 ) =X ( L ) 

26 0 < I .5 ) = 1.0 

CALL MATS ( Q ,CO > N,M , NS I NG ) 
IF (NSING) 72,72,10 

72 CONTINUE 
A=CO (1,1) 

B=CO (2,1) 

C=CO( 3,1) 

D=CO ( A , 1 ) 

E = 0 ♦ 

KANT =2 

GO TO 19 
10 N = 4 

DO 13 1=1, N 
L=K+ 1-1 

Q ( I , 1 ) =X ( L ) **2 
Q( I»2l=X(L)*Y(L) 

Q ( I,3)=Y(L)**2 
Q ( 1,4) =Y ( L ) 

13 Q ( I ,5 ) =1 • 

CALL MATS (Q, CO, N,M, NSING) 
IF (NSING) 73,73,23 

73 CONTINUE 
A=CO (1,1) 

B=CO (2,1) 

C=CO (3,1) 

D=0. 

E=C0 (4,1) 

KANT =3 
GO TO 19 

31 N = 5 

DO 32 1=1, N 
L=K+I— 1 

Q(I,1)=X(L)**2 

G)(I»2)=X<L>*Y(L> 

Q ( I » 3 ) =Y ( L ) **2 
Q ( 1 , 4 ) =X ( L > 

Q( I,5)=Y(L> 

32 Q(I,6)=1.0 

CALL MATS (Q, CO, N,M, NSING) 
IF (NSING) 74,74,25 

74 CONTINUE 
A=C0 (1,1) 

B=CO (2,1) 
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C = CO< 3* 1 ) 

D=C0(4,1 ) 

E=C0 l 5 » 1 ) 

K.ANT = 4 
19 CONTINUE 

BB=B*DVAL ( JGO+O 

CC=C*DVAL ( JGO)**2+E*DVAL( JGOJ-l. 

DSCRIM=BB**2-4**A*CC 

IF CDSCRIM) 7i8>8 

8 FOR T= SORT ( DSCRI M ) 

VI = ( -BB+FORT ) / ( 2 # 0*A ) 
V2=t-BB-F0RT ) / ( 2# 0*A ) 

GO TO 11 
7 CONTINUE 

GO TO (9*10,23*25) » KANT 

9 VARY=-.5*D/A 
GO TO 51 

53 VARY=VARY+FINC*R1< J60 ) 

51 CONTINUE 

VARY=VARY/R1 ( JGO) 

VAL=VAL/R2( JGO) 

56 CONTINUE 

ICNT ( JGO ) - I COUNT 

Q2=VARY 

Q6=TE5T 

RETURN 

54 TE5T=0.0 
GO TO 56 
END 
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$1 BFTC FMXMN 

CFMXMN SUBROUTIME FMXMN 

SUBROUTINE FMXMN ( VAL* VARY >FMNMX ,FINC * S TP INC , TOL ,TEST * J ) 
DIMENSION X<15 > >Y(15) »Q(5,6) >CO<5>2) ,ICNT<3) ,wl<3> ,W2(3> 
DIMENSION FKTR ( 3 ) ,R1(3)»R2(3) 

I COUNT = I CNT < J ) 

Z1=W1 i J ) 

Z2=W2( J) 

K 1 = 5* J— 4 
K2=K1+1 
K3=K2+ 1 
K4=K3+1 
K5=K4+1 

IF (TEST) 5*3*5 

3 SC1=1.0 
Z1=0.0 
Z2 = 0 *0 
TEST -1*0 

I COUNT =— 1 
M=1 

R1{J)=ABS(1 #/VARY ) 

R2( J)=ABS< l./VAL) 

X ( K1 ) =VARY*R1 ( J ) 

Y(K1 )=VAL*R2( J) 

VARY =X ( JC1 ) 

VAL = Y< K1 ) 

GO TO 52 

5 SC1= FK.TR ( J ) 

VARY=VARY*R 1 ( J ) 

VAL= VAL*R2 ( J) 

IF ( ICOUNT-3 } 4*6,6 

4 I COUNT = I COUNT + 1 

6 IF { Z2 ) 35*7,7 

7 K=5*J 
K0=K-4 

1 X(K)=X(K-1) 

Y(K)=Y(K-1) 

K=K~1 

I F ( K-KO) 2,2*1 

2 X ( K1 ) =VAR Y 
Y ( K1 ) =VAL 

IF (Zl) 9*9*8 

8 SC 1-SC1 /2 * 0 
21=0.0 

9 IF (FMNMX*( Y(K1)-Y(K2) ) >10*10*50 

10 Y ( K1 ) -Y < K2 ) 

Y(K2 )=VAL 
X(K1 )=X(K2> 

X(K2 )=VARY 

I F { Z2 ) 11*11,12 

11 SC1=-2.0*SC1 
21=1.0 
Z2=1.0 

GO TO 52 

12 Z2=~l 

GO TO 22 
35 CONTINUE 
DO 55 1=1,4 
L=K1+I-1 
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CO ( I * 1 > =X ( L ) 

35 CO< I » 2 ) =Y ( L ) 

IF ( ( X ( K1 ) -X( K2 > ) / ( y ( r 1) -VARY ) ) 14,14,17 

14 IF C FriNMX* ( VdL~Y ( <1 ) ) ) 16,16,16 

15 CO ( 2 » 1 ) =X ( <1 ) 

CO { 2 ,2 ) =Y ( < 1 ) 

C0(4,1)=X(K2) 

CO ( 4 * 2 ) =Y ( K2 ) 

GO TO 19 

16 CO [ 4 » 1 ) =X ( K3 ) 

C0(4,2)=Y<<3) 

IF ( FMNMX* ( VAL— Y ( k2 ) ) ) 56,56,57 

56 CO (3,1) =VAR Y 
CO 13,2) =VAL 
GO TO 21 

57 CO (2,1) =VARY 
CO ( 2 » 2 ) =VAL 
CO ( 3 » 1 ) =X ( K2 ) 

COO, 2) =Y(K2) 

GO TO 21 

17 I F ( F MNMX* ( VAL-Y ( K1 ) ) ) 20,20,18 

18 CO ( 3 > 1 ) =X ( K2 ) 

CO ( 3 » 2 ) =Y ( K2 ) 

CO ( 2 » 1 ) =X ( K1 ) 

CO ( 2 ,2 ) =Y ( K1 ) 

CO ( 4 » 1 ) =X ( K3 ) 

CO ( 4,2 ) =Y ( K3 ) 

19 COtl,l)=VARY 
C0( 1 *2 ) =VAL 
60 TO 21 

20 CO (2,1) =VAR Y 
C0( 2 ,2 ) =VAL 
CO ( 4 , 1 ) =X t K2 ) 

CO ( 4 ,2 ) =Y ( K2 ) 

21 X(K5)=X(K4) 

Y ( K5 ) -Y ( K4 ) 

DO 51 1=1,4 
K I =K 1+1 

X t K I “1 ) =C0 ( 1,1) 

Y(KI-1)=C0( 1,2) 

51 CONTINUE 

22 GO TO (23,25,31 ) ,ICOUNT 

23 N = 3 
KANT =4 

DO 24 1=1,3 
L =K 1 + 1 - 1 
Q< I » 1 >=X( L )**2 
Q ( I , 2 ) =X ( L ) 

0 ( I , 3 ) = Y ( L ) 

24 Q( I,4)=1.0 

CALL MATS (0,CC,N,i v 1*\SIi\G) 

IF (NSING) 70,70,58 
70 CONTINUE 
A-CO ( 1 * 1 ) 

8 = 0. 

C = 0. 

D=CO (2,1) 

E=CO (3,1) 

GO TO 33 
N =4 


25 
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KANT = 1 
DO 26 1 = 1*4 
L=K1+I-1 
Q( I *1 ) -X ( L ) 

Q ( I *2>=X(L )*Y( L ) 

Q ( I , 3 } =Y ( L ) 

Q ( I>4)=X( L) 

26 Q { I »5)=1*0 

CALL MATS ( Q* CO ,N ,M * NSI NG ) 
IF fNSING) 73*73,29 
73 CONTINUE 
A=CO ( 1 * 1 ) 

B=CO (2*1) 

C=CO(3*l) 

D=CO( A * 1 ) 

E = 0. 

GO TO 33 

29 N = 4 
KANT = 2 

DO 30 1=1*4 
L=K1+ I -I 
0( I *1 )=X(L)**2 
Q ( I * 2 ) =X ( L ) *Y ( L ) 

Q< 1 , 3 ) = Y { L ) **2 
Q{ I *4 ) =Y ( L ) 

30 Q { I * 5 ) = 1* 0 

CALL MATS ( Q * CO *N , M ♦ NSI NO ) 
IF I NSI NG ) 71*71,23 

71 CONTINUE 
A=CO( 1,1) 

B=CO (2,1) 

C=CO (3,1) 

D=0.0 

E=CO (4*1) 

GO TO 33 

31 N=5 
KANT = 3 

DO 32 1=1*5 
L=K1+ I -I 
Q( 1*1 )=XCL>**2 
Q ( I *2)=X(L)*YCL) 

Q( 1*3 )=Y£L)**2 
Q t I * 4 ) =X £ L ) 

G( I *5 ) =Y { L ) 

32 Q ( I *6 ) = 1 « 0 

CALL MATS ( 0, CO ,N ,M * NSI NG ) 
IF C NSI NG ) 72*72,25 

72 CONTINUE 
A= CO ( 1*1) 

B=CO< 2,1) 

C = CO f 3 * 1 ) 

D= CO (4*1) 

E = C0C5* 1 ) 

33 CONTINUE 
AA=C-.2^*S**2/A 
BB=E-.5*B*D/A 
CC=-.25*D**2/A-1*0 
IF t AA ) 27*37*27 

27 QUAN=6B**2-4.*AA*CC 
IF (QUAN) 34,38,38 
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34 GO TO <29,23,25,58 ) »KANT 
. 37 YPRED=-CC/BB 
GO TO 39 

38 YPRED=.5*<-BB+SQRT<QUAN> ) /AA 
KAN=1 

GO TO 39 

28 YPRED=.5*(— B8-SQRT (QUAN J J / AA 
KAN=2 

39 VAREE=-.5*<B*YPRED+D)/A 
GO TO 60 

58 TERM1=( <Y(K3)-Y(K1 ) ) / < X ( K3 )-X l K1 ) ) - ( Y t K. 1 ) -Y C K 2 > ) / < X < K1 )-X C K2 ) ))/ 

1 (X(K3)-X(K2) J 

TERM 2= ( Y ( K1 ) *-Y ( K2 ) ) / ( X ( K1 )-X ( K2 ) 3 
VAREE=. 5* ( X ( K2 ) +X ( K1 ) -TERM2/TERM1 ) 

YPRED=TERM1*VAREE**2+(TERM2-(X<K2)+X(K1>)*TERM1)*VAREE+Y(K2)“X<K2) 
1 *TERM2+X ( K 1 )#X ( K2 ) *TERM1 
60 CONTINUE 

WRITE (6*100)A»B,C*D»E»BB*CC> 

100 FORMAT (2H0A,E15.7»3H B,E15.7,3H C»E15.7,3H 0.E15.7.3H E,E15 

1.7.4H BB,E15.7,4H CC,E15.7) 

VARI=VAR£E*57.295775/R1 (J) 

Y1PRED=YPRED/R2 < J) 

WRITE 16,200) Y1PRED, VARI 
200 FORMAT (2E18.7) 

I F ( X ( K1 )-X ( K2 > ) 41,41,42 

41 SIG =1.0 
GO TO 43 

42 SIG =-1.0 

43 IF (SIG*( VAREE-X1K2) ) 145,44,44 

45 I F { X ( K.1 )-X ( K3 ) ) 46,46,47 

46 SIG =1.0 
GO TO 48 

47 SIG =-1.0 

48 IF (SIG*( VAREE-X1K3) ) >59,44,44 

44 60 TO (28, 34), KAN 

59 IF ( ABS<VAL-YPRED)-T0L*R2(J) ) 54,54,13 
13 VARY=VAREE 

GO TO 53 
50 Z2=1.0 

SC1=STPINC*SC1 

52 VARY=VARY+SC1*FINC*R1C J) 

53 ICNT 1 J)=ICOUNT 
FKTR ( J)=SC1 
Wl( J )=Z1 

W21 J)=Z2 
VARY=VARY/R1(J) 

VAL=VAL/R2 l J ) 

RETURN 

54 TEST=0.0 
GO TO 53 
END 



